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INTRODUCTION 

Research during the period June 1974 through May 1975 has resulted in 
the following Fortran language computer codes: (1) more efficient two- and 

three-dimensional central force potential solvers; (2) a three-dimensional 
simulator of an isolated galaxy which incorporates the aforementioned poten- 
tial solver; (3) a two-dimensional particle-in-cell simulator of the Jeans 
instability in an infinite self-gravitating compressible gas; and (4) a 
two-dimensional particle-in-cell simulator of a rotating self-gravitating 
compressible gaseous system of which rectangular coordinate and superior 
polar coordinate versions were written. 

POTENTIAL SOLVERS 

The two- and three-dimensional potential solvers listed in Appendices 
A and B respectively decrease computing cost as compared to the previous 
versions by 50 and 75 percent respectively. Their methods of operation 
were detailed in the Semi-Annual Status Report for the period June 1974 
through November 1974. The three-dimensional potential solver (pages B-4 
through B-17 of this annual report) is virtually self-explanatory due to 
extensive use of comment cards. 

THREE-DIMENSIONAL ISOLATED 
GALAXY SIMULATOR 

Through multiple overlays the three-dimensional potential solver has 
been incorporated into an efficient n-body dynamic code, which provides for 
plotting of various views of the particle distribution in position and velocity 
space as well as printed diagnostics. A set of balanced initial conditions have 
been tested for a thin disk of particles. In creating the initial conditions. 
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the particles are first randomly distributed according to an axi-symmetric 

2 1/2 

radial density distribution which varies as [1 - (r/r Q ) ] , where r 

and r are the radius and maximum radius respectively, and according to 
o 

a Gaussian axial density distribution. After solving for the potential, 
the particles are assigned a Maxwellian velocity distribution after deter- 
mining the component standard deviations. That of the axial velocity is 
determined by balancing the axial components of the gradients of the potential 
and the pressure. The standard deviation of the radial and azimuthal velocity 
distributions are determined by dividing the galaxy into layers of constant 
axial position and by applying to each layer the Toomre stability criterion 
for infinitesimally thin disks. An average azimuthal velocity is calculated 
for each radius and height by reducing the square of the angular velocity of 
the balanced cold disk by terms involving (a] the radial pressure gradient, 

(b) the difference in the squares of the standard deviations of the radial and 
azimuthal velocities, and (c) the axial derivative of the expectation value of 
the product of the radial and axial velocities. 

This code is listed in Appendix B; position and velocity space plots of 
a set of initial conditions and the first few cycles of a run are presented 
in Appendix C. 

Longer runs on larger meshes (up to 64 x 64 x 16) with initial conditions 
computed with this method and probably other methods will be made in the near 
future . 


*Toomre, Alar: On the Gravitational Stability of a Disk of Stars. 

Astrophys. J., vol. 139, no. 4, May 15, 1964, pp. 1217-1238. 
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TWO -DIMENS TONAL GASEOUS JEANS 
INSTABILITY SIMULATOR 


The particle-in-cell method attempts to combine the advantages of the 
Lagrangian and Eulerian approaches to non-steady compressible gas simulation. 

The basic method as described in Ths ParticIe-in-Cell Method for the Calcula- 
tion of the Dynamics of Compressible Fluids by A. A. Amsden, Los Alamos 
Scientific Laboratory Report LA-3466 of February 1966, was modified to 
include self-gravitation and periodic boundary conditions. Self-gravitation 
was implemented by adding gravitational terms to the cell calculation of 
new velocity and specific internal energy. Periodic boundary conditions were 
implemented by (a) requiring that particles leaving one boundary carry their 
mass, energy, and momenta into the opposite boundary and (b) finite differencing 
pressure and potential near the boundary in such a way that a value required 
from a cell just outside the boundary is obtained from a cell just inside the 
opposite boundary. 

This code is listed in Appendix D and plots from two full length runs 
(150 cycles) made on a small mesh (32 x 32 cells) are presented in Appendix E. 
These runs verify the analytic prediction that an increase in initial tempera- 
ture decreases both the number of instabilities per unit area and their rate 
of evolution. That the code conserves total energy is evident from the 
plots of Appendix E while printed diagnostics (not shown) verify conservation 
of linear momentum. 

When point masses of a collisionless gravitational n-body code were given 
an initial position distribution identical to that given the gas "particles" in 
this particle-in-cell fluid code and an initial velocity dispersion corresponding 


to the initial specific internal energy of the cells in this fluid code, an 
almost identical pattern of instabilities resulted. The dynamics of the 
instabilities occuring in these two codes will be compared in the near 
future. Also planned is the combination of these two codes for use in 
studying stellar-gaseous gravitational two-stream instabilities. 

POLAR COORDINATE SIMULATOR FOR 
A ROTATING GAS 


The status report for June 1974 - November 1974 described the modifica- 
tion of a two- dimensional particle-in-cell code wr.^h square cells and linear 
dynamics to a code with square cells but with angular and radial momentum 
transfer during the particle movement phase. Although this modification 
conserved angular momentum and greatly reduced artificial heating, extensive 
testing demonstrated that the combination of square cells and polar dynamics 
produced a large assymetric radial outward acceleration. A "quick fix" of 
rotating the meshes (+) or (— ) 45 on alternate cycles served only to symmetrise 
the artificial radial acceleration. 

Recently a code was written which combines a rectangular mesh potential 


solver with polar meshes for the cell values of radial and angular momentum, 
specific internal energy, mass, x and y coordinates of the cell center of 
mass, and the average radius squared for the center and end of each time step. 


At the center of mass of each polar cell a combination finite differencing and 
bilinear interpolation of the rectangular potential mesh yields the x and 
y components of gravitational field which are then resolved into radial and 
gravitational components . The particle movement phase depends on the princi- 
ple that for a collection of particles moving with a uniform angular velocity, 
e.g., a cell, both the rotational kinetic energy and the angular momentum are 


proportional to the average of the radius squared (radius of gyration squared) times; s 




a power of the angular velocity. During the particle movement phase, the new 
» 3s density is built up simultaneously in the polar mass density mesh and in 
the rectangular mesh which was used to store the potential during the cell 
calculations. The possibility of high random cell velocities near the center 
of the gas was provided for by making the radial and azimuthal dimensions of 
all cells roughly equal to one unit; this was accomplished by increasing the 
number of cells per mass ring from four for the first ring to 128 for the 
31st. ring. Since the basic particle-in-cell method is a leap frog method, 
it is critically important that all positions, forces and accelerations be 
computed at the center of the time step and that all velocities, momenta, 
kinetic energies and the associated radii of gyration be computed at the end 
of the time step; a verification of this code is that identical results were 
produced by alternate procedures for leap frogging velocity and radius of 
gyration one half time step ahead of particle position. 

The above polar coordinate code is listed in Appendix F. Appendix G 
presents plots produced by runs of the polar and earlier rectangular codes, 
describes their identical initial conditions, and demonstrates the superiority 
of the polar code by a brief analysis. The polar code suffers from a very 
small radial outward acceleration and this is being investigated with a view 
toward elimination. 

To conserve computer resources, the system state at the end of a run 
may be stored on magnetic tape and the run may later be continued if desired. 
If the time scaled velocity of any cell exceeds the dimension of that cell, 
the run is terminated with a complete set of plots and a long printout of 
cell quantities and a continuation tape is generated for possible later 
use with a smaller time step. 
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Plans for the use of this code include investigation of the evolution 
of (a) a purely gaseous system [listing of Appendix F with one parameter 
change (b) a gaseous system acting under the influence of an analytically 
computed con: :.ant stellar central force component (Appendices F and G) and 
(c) a gaseous system interacting through the gravitational potential with 
an existing collisionless stellar n-body code. The code may be modified 
to allow star formation and/or permit internal energy loss by electro- 


magnetic radiation. 


APPENDIX A 


Listing of the Two-Dimensional Potential Solver of Increased Efficiency. 

This subroutine decreases cost by 50% by replacing central memory 
storage with disk file storage. To minimize file buffer size and peripheral 
process time it utilizes the input/output routine MEMDISK and associated 
assembly language subroutine MDFUNC which were written by R. Bulle of the 
Analysis and Computations Division of the Langley Research Center. 


Subroutine Name 

Page No. 

GETPHI 

A- 2 

MEMDISK 

A- 9 

MDFUNC 

A-ll 
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SUBROUTINE GETPHI 

COMMON Z < 1025) *Y( 1025) , RHO U2B* 1 28 ) * I2A* lTEST*F ( 8385 ) 
DIMENSION RHO 1 (128*64) , RH02 < 1 28* 64) * HI N2 1(1 28 ) *HJN21 (128) 
DIMENSION IPAR (5) 

EQUIVALENCE (RHO ( 1 . I ) » RHOl ( 1 » 1 ) ) * ( RHO (1 » 65 ) , RH02 ( 1 * 1 ) ) 
IF(ITEST.EQ.O) GO TO 70 

itest=o 

I2B= I2A-1 

N=2*»I2A 

N02=N/2 

N21=N02*1 

M04=N/4 

N34=N02*N04 

NN24=N02»N04 

RNI=1 ./ (N*N) 

DO 5 J=1*N02 
DO 5 1=1* N02 

IF(I.EQ.l.AND.J.EQ.I) GO TO 5 
R2=<I-I)*(I-1)+(J-1)*<J-1) 

R=SQRT(R2) 

RHO(I, J)=RNI/R 
5 CONTINUE 

RHO ( 1 ♦ 1 ) =RHO ( 1 *2 ) 

N02SQ=N02<»N02 
DO 10 J=1,N02 
R2=N02SQ+( J-I)*(J-1) 

R=SQRT(R2) 

HIN21 ( J) =RNI/R 
10 HJN21 (J) =HIN21 (J) 

R2=N02SO*N02SQ 
R=SQRt (R2) 

HI JN2 1 =RNI/R 
CALL GETSET ( 2* I2B ) 

DO 20 J=1 *N02 
DO 15 1=1 *N02 
15 Z ( I ) =RH0 ( I * J) 

Z (N21 ) =HIN21 (J) 

CALL FTRANS (2* I2B) 



20 


HIN21 { J) =Y (N21 ) 

DO 20 1=1, N02 
RHO { I , J ) =Y ( I ) 

DO 25 1=1 »N02 
25 Z (I)=HJN21 (I ) 

Z (N21 ) =HI JN21 
CALL FTRANS(2,I2B) 

HI JN21=Y (N21 ) 

DO 30 1=1, N02 
30 HJN21 ( I ) =Y ( I ) 

DO 40 1=1, N02 
DO 35 J=1 ,N02 
35 Z<J)=RH0U,J> 

Z(N21)=HJN21 (I) 

CALL FTRANS (2, I2B) 

HJN21 ( I ) =Y (N21 ) 

DO 40 J=1,N02 
40 RHO( I , J) =Y { J) 

DO 45 J=1 ,N02 
45 Z(J)=HIN21 (J) 
Z(N21)=HIJN21 
CALL FTRANS ( 2 » I2B ) 
HIJN21=Y(N21) 

DO 50 J=1 , N02 
50 HIN21 (J)=Y{J) 

DO 55 1=1, N02 
DO 55 J=1,I 

M=lMj-l)*N21-U-l)*J/2 
55 F (M)=RHO ( I , J) 

DO 60 J=1,N02 
H=N21 ♦ ( J-l ) fr N21- ( J*"l > *J/2 
60 F (M)=HIN21 { J) 

NFM=N21 +N02*N2l-N02«N21/2 

F <NFM) =HI JN21 

IPAR (1 ) =1 

IPAR (2) =0 

IPAR (3) =NFM 

IPAR (4)=5LTAPE5 

CALL MEHDISK(F(1>,IPAR) 



REWIND 5 
GO TO 200 
70 IPAR ( 1 > =0 
IPAR(3)=NFM 
IPAR(4)=5LTAPE5 
CALL MEMDISK(F(I)*IPAR) 
REWIND 5 
72 IPAR ( 1 ) =1 
IPAR<3)=NN24 
IPAR (4) =5LTAPE3 
CALL MEMDI SK (RH02 (1,1) » IPAR ) 
REWIND 3 

CALL GETSET (3,I2A) 

DO 80 J= 1 * N04 
DO 75 1=1, N02 
Z (I ) -RHOl ( I » J) 

75 Z (N02* I ) =0 • 

CALL FjRANS {3» I2A) 

DO 80 1=1, N02 
RHOl <I,J)=YCI) 

80 RH02 { I * J ) =Y ( N02* I ) 

IPAR ( 1 ) -1 

IPAR (4)=5L TAPE 1 

CALL MEMDISK(RH01 (1,1) , IPAR) 

REWIND 1 

IPAR {4 ) =5LTAPE2 

CALL MEMDI SK <RH02 (1,1), IPAR) 

REWIND 2 

IPAR ( 1 ) =0 

IPAR(4)=5LTAPE3 

CALL MEMDISK (RHOl (1,1), I PAR) 

REWIND 3 

DO 90 J=1 , N04 

DO 85 1=1, N02 

Z<I)=RH01 (I,J) 

85 Z (N02+ I ) =0 • 

CALL FTRANS (3, I2A) 


DO 90 1=1, N02 
RHOl ( I * Jl =Y ( I ) 

90 RH02 ( I, J>=Y (NO 2+1) 

IP AR ( 1 ) = 1 

I PAR (A) =5LT AREA 

CALL MEHDISK(RH02U*15 ,IPAR) 

REWIND A 

IPAR { 1 } =0 

IPARU)=5LTAPE1 

CALL MEMDISK(RH02{1,1) ,IPAR) 

REWIND I 

DO 115 1=1, N02 

DO 95 J=1 »NO& 

2 ( J) =RH02 ( I , J1 
Z<N04,J)=RHOl <I,J> 

Z (N02+ J) =0 • 

95 Z I N34+ J) =0 • 

CALL GETSEt (3, I2A) 

CALL FTRANS (3, I2A) 

DO 100 J=2,N02 
IFU.LT.JISO TO 96 
M=I+ ( J-l )*N21- ( J“1 ) °J/2 
GO TO 97 

96 M=JMI-l)*N21-U“l>*I/2 

97 Z (J)=Y(J)*F(M) 

100 Z (N02+J) =Y (N02+J) *F (M) 
ZC1>=Y(1)*F(I) 

M=N21 +< 1-1) »N21- <1-11*1/2 
Z(N21)=Y<N21)*F<M> 

CALL GETSET (A, I2A) 

CALL FTRANS (A, I2A) 

DO 115 J=l,NOA 

RHOl ( I , J ) =Y ( J ) 

115 RHQ2 ( I , J) =Y ( NOA+J) 

IPAR ( 1 ) =1 

IPAR(41=5LTAPE1 

CALL HEMDISK(RH01 (1,1) , IPAR) 



REWIND 1 

I PAR (4-) =5LT APE 3 

CALL MEMDI5K (RH02 ( 1 * 1 ) * IPAR ) 

REWIND 3 

IPAR ( 1 ) =0 

IPAR (A) -5LTAPE2 

CALL HEM0ISK<PH01 ( 1 » 1 > • IPAR ) 

REWIND 2 

IPARU)=5LTAPE4 

CALL MEMDISK (RH02(I*1> tIPAR) 

REWIND 4 

DO 130 I“1 *N02 

DO 120 J=1*N04 

Z(J)=RH01 (I»J) 

Z (N04+ J ) =RH02 ( I * J) 

Z (N02+ J) -0 • 

120 Z(N34+J>=0. 

CALL GETSET ( 3* I2A) 

CALL FTRANS ( 3* I2A) 
IF(I.EQ.l)GO TO 125 
DO 123 J=2»N02 
IF(I.LT.J)GO TO 121 
N=I+ <J-l>*N21-<J~l)*J/2 
GO TO 122 

121 (I-l)*N2l-(I-l)*I/2 

122 Z(J)=Y<J>*F*M> 

123 Z (N02 + J > = Y { N02+J) *F (M ) 
ZU)=Y(1)*F(I> 

H-N2 l+(I-l)*N21-(I-l)#I/2 
Z(N21)=Y(N21)*F(M) 

GO TO 128 

125 DO 127 J=2*N02 

M=N21+ { J~1 ) *N21- ( J~1 ) *J/2 
Z(J)=Y(J)*F(M) 

127 Z(N02*J)=Y(N02+J)*F(M) 
ZU)=Y(l>*F<N2l> 
Z<N21)=Y{N21)*F(NFM) 



128 CALL GETSET(A,I2A) 

CALL FTRANS (A* 12A) 

DO 130 J=l*NOA 
RHO 1 < I ♦ J ) =Y ( J ) 

130 RH02 ( I * J) =Y (NOA+J) 

IPAR { 1 ) — 1 

IPAR ( A ) =5LTAPE2 

CALL MEMDI Sk IRHO l < 1 * 1 ) » IPAR) 

REWIND 2 

IPAR ( 1 ) =0 

IPAR (A) =5LTAPE3 

CALL MEMDISKIRH01 (1*1) , IPAR) 

REWIND 3 

DO 1 AO J = 1 *N0A 

DO 135 1=1 *N02 

Z ( I ) =RH01 ( I * J) 

135 Z (N02* I ) =RH02 ( I * J) 

CALL FTRANS (A* I2A) 

DO 1 AO 1=1 *N02 

1A0 RHOl ( I * J) =Y < I ) 

IPAR ( 1 ) =1 

IPAR (A)=5LTAPE3 

CALL HEMDISK(RH0I(1*1)*IPAR) 

REWIND 3 

IPAR ( 1 ) =0 

IPAR (A) =5LTAPE1 

CALL MEMDISK(RH01 (1*1) *IPAR) 

REWIND 1 

IPAR (A) =5LTAPE2 

CALL MEM DISK (RH02 (1*1) *IPAR) 

REWIND 2 

DO 150 J=1*N0A 

DO 1A5 1=1 *N02 

Z ( I ) =RH01 ( I * J) 

1A5 Z (NQ2+I ) =RH02 ( I * J) 

CALL FTRANS (A, I2A) 



DO 150 1-1 t N02 
150 RH01 ( I ♦ J) = Y ( I ) 

IPAR ( 1 ) =0 
IPAR(4)=5LTAPE3 
CALL MEMDISK<RH02(1,1),IPAR) 
REWIND 3 
200 RETURN 
END 



oooooooonoo 


SUBROUTINE MEM0ISk<I8UFF,IPARAMS) 

DIMENSION JPAR (5) 

equivalence 
I (JPAR { I ) ♦ IFUNCT )* 

2 (JPAR (2) t IMODE >» 

3(JPAR{3) , LENGTH >♦ 

4 (JPAR (A) * IAGTN )« 

5( JPAR (5) * I STATUS) 

DIMENSION IPARAM5 (5) 

DIMENSION IBUFF ( 1 3 0 > 

«■**«■«■ •ft********** «■•»*****■» 

IBUFF - BUFFER 
IPARAMS - PARAMETER LIST 

1 IFUNCT = FUNCTION 0=READ ti-WRlTE 

2 IMODE = MODE 0=BINARY « l=COD£D 

3 LENGTH = LENGTH OF MESSAGE 

5 iStatus s error status 

0 NO ERROR 

1 invalid mode 

2 invalid function 

IFUNCT =IPARAMS(1) 

IMODE =IPARAMS(2) 

LENGTH -IPARAMS (3) 

IAGTN = IPARAMS (4) 

ISTATUS=0 
IF ( IMODE) 44*45*43 

43 IF ( IMODE-1 >45»45»44 

44 IStATUS=1 
GO TO 999 

45 CONTINUE 
IF ( IFUNCT) 9 *7*B 

7 IFUN=1 OB 

IF(IMODE.EQ.O) IFUN=12B 
GO TO 11 

8 IF( IFUNCT-1 )9*10»9 

9 ISTATUS=2 
GO TO 999 

10 IFUN=34B 

IF(IM0DE.EQ«0> IFUN=36B 
■>- 11 CONTINUE 
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IOENT 

MDFUNC 


entry 

MDFUNC 


ENTRY 

mdsetf 

RAPLUSl 

MACRO 


+ 

SA1 

RA + 1 


NZ 

endm 

Xl»* 

RA 

EQU 

0 

FN 

BS5 

0 

FET 

DATA 

0 

fst 

data 

0 

IN 

data 

0 

OUT 

data 

0 

LIM 

data 

0 

CIOCALL 

vfd 

18/3LCI0* I/O* 1/1 *22/0* 18/FET 


vfd 

42/0HMDFUNC*l8/3 

hofunc 

* 

DATA 

0 


* CALL MDFUNC(IFUNCT*L*STAT> 


« 


SA1 

B1 

,xi=ifunct 

MXO 

42 

,X0=42BITS 

SA2 

FET 

•X2=FET 

BX3 

-X0*X1 

• X3=l 8R( IFUNCT ) 

BX4 

X0*X2 

.X4=:42U(FET) 

SB7 

1 

• B7=l 

BX6 

X4 + X3 

.X6=42U(FET> ♦ ISR(IFUNCT) 

SA6 

A2 

.FFj=:FET + IFUNCt 

sai 

CIOCALL 

.Xl=CIO 

BX6 

XI 

• X6=CI 0’ 

RAPLUSl 

•lSRA+1 CLEAR 

SA6 

B7 

.CALL CIO 

XJ 


•exchange jump to monitor 

♦ Sai 

A2 

.xi=fet 

SA2 

IN 

.X2=IN 

BX6 

-X0*X1 

.x6=code and status 

SA3 

A2*B7 

.X3=0Ut 

SA6 

B3 

.STAT = CODE AND ST ATUS 

1X7 

X2-X3 

.X7=IN-0UT 

PL 

X7*FFT1 

•IF(IN.GE.OUt>GO TO FFTl 


r> 

\ 


■ 



! SAl A2-B7 .Xl=FIRST 

f SA4 A3 + B7 .X^LIHIT 

I 1X7 X7-X1 .X7=IN-FIRST-0UT 

1X7 X7 + X4 .X7=IN-FIRST*LIHIT-0UT 

FFTl Sa7 B2 .L = IN-FIRSt<-LIMIT-OUT 

.1 EQ MDFUNC 

f VFD 42/OHMDSETF * 18/2 

mdsetf DATA 0 

;$ ’ •» 

* CALL MDSETF ( I BUFF »JPAR AMS) 


>'■ 


SB6 

3 

1. 


SB4 

1 

;?• 


Sax 

B2 + B6 

?. 

1 


SB5 

B4 + B4 


SB3 

FET 

1 

ft 


BX6 

XI 

k 


SB7 

B4 + B6 

§ 


SA6 

B3 



SAl 

B2 



SX6 

B1 

4 


ZR 

Xl *READ 

!• 

i 

WRITE 

SA5 

B2 + B5 

•Vi 


1X7 

X6*X5 

T 


SA6 

B3 + B6 

I 


SA7 

B3 + B5 

5:-' 

1 


EO 

BOTH 

r ,5 

READ 

SA6 

B3 + B6 

1 


SA6 

B3 + B5 

> 

fl 

BOTH 

SA5 

B2 + B5 

1 ' 


5X4 

B4 

'V. 


1X7 

X6 + X5 

|. 


SA6 

B3 + B4 

% c 


1X7 

X7 + X4 



SA7 

B3*B7 



EQ 

END 

mdsetf 


' H 

iti- 

r'iS— 

I3> 

* 

IP 


.86 = 3 

• B4=l 

•X 1 =AGT NAME 

• B5=2 
.B3=FET 
•X6=AGT NAME 
.87 = 4 

.FET=AGT NAME 

.xi=ifunct 

.X6=IBUFF 

.IF(IFUNCT.EO.O)GO to read 
.X5=LENGTH 
•X7= I BUFF + LENGTH 
,0UT=IBUFF 
.IN =IBUFF-vLENGTH 
.GO TO BOTH 
.OUT=IBUFF 
.IN =IBUFF 
• X5=LENGTH 
• X4= 1 

.X7=IBUFF*LENGTH 

.1St=IBUFF 

.X7=IBUFF+LENGTH+I 

.LlM=IBUFF*LENGTH+l 

.RETURN 1 



APPENDIX B 


Listing of the Three-Dimensional Isolated Galaxy Dynamic Simulator. 

The initial conditions are for a thin balanced disk the velocity 
dispersion of which satisfies the Toomre local stability criterion. 

Overlay programs GETH and GETPHI constitute an improved three- 
dimensional potential solver which decreases cost by 75% by replacing 
central memory storage with disk file storage. 

To make the listing easier to understand, use of the input/ output ■ 
routine listed in Appendix A has been replaced in this listing by the less 
efficient Fortran binary READ and WRITE statements and their accompanying 
larger disk file buffers. 


Overlay No. 

Program/Subroutine Name 

Page No. 

(0,0) 

PROGRAM STARS 3D 

B-2 

(1,0) 

PROGRAM GETH 

B-4 

(2,0) 

PROGRAM GETPHI 

B-6 


SUBROUTINE ANLX 

B-9 


SUBROUTINE ANLSYN 

B-ll 


SUBROUTINE SYNX 

B-16 

(3,0) 

PROGRAM INITIAL 

B-18 

(4,0) 

PROGRAM STARS 

B-33 


*Toomre, Alar: On the Gravitational Stability of a Disk of Stars. 

Astrophys. J., vol. 139, no. 4, May 15, 1964, pp. 1217-1238. 


OVERLAY (IFILE»0»0) 

PROGRAM SiARS3D( OUTPUT *TAPE1» TAPE2»TAPE3»TAP£4.TAPE5,TAPE6»TAPE7* 

1 TAPES. TAPER, TAPE10»TAPE11*TAPE12) „ irv , 

COMMON/ ALLC0M/N»N02»N21 .N04.N41 »N34,NH»NH02,NH21 »NCH,NRH0»NHH» 


1 I2A,I2B,I3A,I3B,NH04 

COMMON/HN21COM/HN21 (33,9) 

COMMQN/INITCOM/MASSD (16,4) 

COMMON/ ADVC0M/NBR»N8S,R I » XM»DT , DTE2 *N04Ml ,N04M2, VXYMAX » VZMAXI ♦ 

1 CY,CYY,NPLOT,NPRINT» IN (2) ,XMIN»XMAX,YMIN,YMAX,ZMIN, 

2 ZMAX,RMIN,RMAX,VRMIN,VRMAX,VTMIN,VTMAX*VZMIN»VZMAX,XPP (3> , 

3 Yp,ZP, YPP (3) ,RPP (3) ,VTP, VRP, VZP, ITAPX,PI,MASK1 ,S2, JT, J5, NMK, 

4 NBS3,CXY,CZ,VZAV1,DR,ItES T ,JTFTLE,JSFILE,DDD,N04P1 

integer cy.cyy 

REAL MA5SD 

C NAME OVERLAY FILES A 

C NAME THE FILE FOR THE GETH AND INITIAL OVERLAYS 

IFILE=5LIFILE 

C NAME GETPHI OVERLAY FILE 

GPFILE=6LGPFILE 
C NAME STARS OVERLAY FILE 

SFILE=5LSFILE 
C SET DIMENSIONS 


I2A=6 


I3A-4 
I2B= I2A-1 
I3B=I3A-1 


N=2»*I2A 

N02=N/2 

N21-N02+1 

N04=N/4 

N04P 1 =N04+ 1 

N41=N04+ I 

N34=N02+N04 

NH=2*»I3A 

NH02=NH/2 

NH2-1=NH02 + 1 

NH04=NH/4 

NCH=N02*N02*NH02 

NRH0=N04*N04»NH02 

NHH=N04*N04«NH21 

C CALL GETH OVERLAY IN ORDER TO COMPUTE THE TRANSFORMED GREENS FUNCTION AND 

c store chunks of it on disk file 9 in the order appropriate for later use in 


t 

i 


-% 



Li?; 


C THE GETPHI OVERLAY 

CALL OVERLAY (IFILE*i*0»6HRECALL) 

C CALL I TH£ T INITIAL OVERLAY IN ORDER TO SET INITIAL VALUES, GENERATE INITIAL s TARi 
C POSITIONS, and compute INITIAL DENSITY. 


CALL OVERLAY ( IFILEf3*0,6HRECALU 
C CALL GETPHI OVERLAY IN ORDER TO COMPUTE THE 


INITIAL POTENTIAL: FROM THE INITIAL! 


C density. 

CALL OVERLAY (GPFILE,2,0»6HRECALL) 

itest=i 


CYY=0 

C CALL INITIAL OVERLAY IN ORDER TO ASSIGN INITIAL VELOCITIES TO STARS, 

c time Step the initial star positions, compute a new initial, density, 
C AN INITIAL SET OF PLOTS. 

CALL OVERLAY (IFILE»3»0,6HRECALL) 


HALF 

AND MAKE 


CYY = 1 

10 PRINT 12, CYY 

12 FORMAT I7H CYCLE=*I8> 

CALL SECOND < T 1 > 

C CALL GETPHI OVERLAY IN ORDER TO COMPUTE POTENTIAL! 

CALL OVERLAY (GPFILE* 2, 0*6HRECALL) 

CALL SECOND (T2) 

TGETPHI=T2-Tl 
PRINT 27* TGETPHI 

27 FORMAT ( 12H FIELD TIME=, E16.8) 

c call stars overlay in order to compute new star positions and velocities and 

c A NEW DENSITY MESH. 

CALL OVERLAY (SFILE»4fO,6HRECALL) 


42 

C IF 
C GO 


CALL SECOND (T3) 


TSTARS.-T3-T2 
PRINT 42, TSTARS 


FORMAT ( 1 9H Star ADVANCE TlME=*E16.8) 

It IS NOT the LAST CYCLE (time STEP) CY INCREMENT THE CYCLE NUMBER CYY AND 

to statement 10 . 

IF (CYY.GE.CY) GO TO 50 


45 


CYY *= CYY + 1 
GO TO 10 
STOP 
END 


50 


OVERLAY (IFILE*1*0) 

c this overlay performs a cosine analysis of the three-dimensional sreens 
C FUNCTION array. IT THEN WRITES CHUNKS Or THIS ARRAY ON DISK FILE 9 IN THE 
C ORDER IN WHICH THEY WILL BE READ INTO THE HH ARRAY DURING THE GETPHI 
C OVERLAY. VALUES FOR I=N/2+l AND J=N/2+l ARE TRANSFERRED TO THE HN'21 ARRAY 
C WHICH IS IN COMMON WITH THE GETPHI OVERLAY. 

COMM0N/ALLCOM/N»NO2*N21 » N04t N41 * N34*NH*NH02* NH21 *NCH* NRHO* NHH* 

1 I2A» I2B* I3A* I3B*NH04 
COMMON/HN21COM/HN21 (33*9) 

COMMON Z < 1025) * Y(1025) 

DIMENSION H ( 33 » 33*9) 

RNI-1 •/ (N^N^NH) 

DO 1 K-l *NH2l 

DO 1 J=1*N21 
DO 1 1=1 »N21 

RI=(K-1)*(K-1) 

IFtRI.LT.l.) RI-1. 

H(I*J,K)=RNI/SQRT(RI) 

1 CONTINUE 

CALL GETSET (2* I2B) 

DO 2 K=1 »NH2l 
DO 2 J=1*N21 
DO 3 1=1 *N21 
3 Z(D S H(I*J»K) 

CALL FTRANS (2* I2B) 

DO 4 1=1 ?N21 
A H(I*J.K)=Y(I) 

2 CONTINUE 

DO 5 K=I »NH21 
DO S 1=1 *N21 
DO 6 J=1 *N21 

6 Z< J)=H( I* J»K) 

CALL FTRANS (2* I2b) 

DO 7 J=1*N21 

7 H(I*J*K)=Y(J) 

5 CONTINUE 

CALL GETSET ( 2* I3B) 

DO 10 J=1 * N21 
DO 10 1=1 *N2l 
DO 8 K=1 * NH21 


8 Z (K) -H ( I » J»K) 

CALL FTRANS ( 2f I3B) 

DO 9 K=1*NH21 

9 H< I» J,K>=Y(K) 

10 CONTINUE 

N0AP1=N41 

WRITE(9) < < <H<I»J,K> *1 
WRITE (9) ( ( (H ( I * JtK ) * I 
WRITE (9) ( ( ( H ( I * J*K ) * I 

WRITE (9) < ( <H(I*J,K) *1 

WRITE (9> ( ((H(I»JfK) »I 

WR ITE (9) (((H(I*J*K)*I 
REWIND 9 
DO 15 K=1 »NH21 
DO 15 1=1, N21 
15 HN21 (I*K)=H(I»N21*K> 

RETURN 

END 



:1,N0A) ,J=l,NOA) ,K=1,NH21) 

=1,N04) ,J=N0AP1,N02) ,K=I,NH21) 
=1,N04) ,J=1*N04) ,K=1,NH21) 
=N0AP1,N02) * J=1 ,NOA) *K=1 ,NH21 ) 
=NOAP1»N02> , J=NOAP1,N02) ,K=1*NH21> 
=N04P1,N02) , J=l,NOA) ♦K=1*NH21) 


c 

c 

c 

c 

c 

r 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


OVERLAY (GPFILE»2*0) 

PROGRAM GETPHI 

THIS OVERLAY SOLVES FOR THE POTENTIAL 1 MESH (DIMENSIONED N/2 BY N/2 BY NH/2) 
DUE TO A DENSITY MESH (DIMENSIONED N/2 BY N/2 BY NH/2) BY DOING A PERIODIC 
ANALYSIS OF THE DENSITY AND THEN A PERIODIC SYNTHESIS OF THE PRODUCT OF THE 
TRANSFORMED GREENS FUNCTION (DIMENSIONED (N/2+1) BY (N/2*l) BY (NH^ + l)) AND 
THE TRANSFORMED DENSITY, FORMALLY SPEAKING, EACH OF THE TRANSFORMS EXCEPT 
THE COSINE ANALYSIS OF THE GREENS FUNCTION, WHICH IS PERFORMED IN THE GETH 
OVER! AY) REQUIRES AN ARRAY DIMENSIONED N BY N BY NH/2. TO REDUCE CORE 

storage this overlay performs these transforms in chunks by the alignment 
of four smaller arrays named rhoi»rho2,rho3» and rhoa, each of which is 

DIMENSIONED N/A BY N/A BY NH/2. THE CHUNKS OF THE LARGER ARRAY NOT IN CORE 
AT ANY ONE TIME ARE STORED ON DISK FILES l THROUGH 8. THE FOLLOWING ARE TWO 
TOP VIEWS OF THE LARGER ARRAY BOTH OF WHICH DESIGNATE THE CHUNKS AS IROW AND 


JCOLUMN. IROW 1 AND 2 OF JCOLUMn 1 AND 2 CONSTITUTE TW E ACTIVE MESr. IN 
THE DIAGRAM ON THE LEFT THE NUMBERS WITHIN THE CHUNK'- • r JCOLUMN 1 AND 2 
INDICATE THE DISK FILES ON WHICH THOSE CHUNKS ARE STuicD. (NO DISK STORAGE 
IS REQUIRED FOR JCOLUMN 3 OR A.) REFERRING TO THE DIAGRAM ON THE RIGHT* THE 
NUMBERS WITHIN THE CHUNKS ARE THE ORDER IN WHICH CHUNKS OF THE TRANSFORMED 

density are multiplied (Element by element) by the appropriate portion of the 

TRANSFORMED GREENS FUNCTION WHICH HAS BEEN READ FROM DISK FILE 9 INTO ARRAY 
HH (N/A » N/A* NH/2+ 1 ) . (AN EXCEPTION. IS THE SET OF TRANSFORMED GREENS FUNCTION 
BOUNDARY VALUES FOR I=N/2*1 AND J~N/2* 1 WHICH REMAIN AT ALU TIMES IN COMMON 
IN THE ARRAY HN21 (N/2+1 ,NH/2+l ) . ) A PLUS IN A CHUNK INDICATES THAT NEW VALUES 
MUST BE READ INTO ARRAY HH BEFORE THAT CHUNK IS MULTIPLIED BY HH. THIS SYStEMi 
MINIMIZES PERIPHERAL PROCESS TIME BY UTILIZING THE PERIODICITY OF THE 


TRANSFORMED GREENS FUNCTION. 


o o no nonooooooonoo onooooooooo ooooo 


TWO TOP VIEW5 OF LARGER MESH (N BY N BY NH/2) - IROW 1 AMD 2 OF JCOLUMN 
1 AND 2 CONSTITUTE THE ACTIVE MESH (N/2 BY N/2 BY NH/2). THE 

directions are x and i - down: on page, y and j - to right on page, 

Z AND K - OUT of PAGE. 


JCOLUMN 
12 3 4 


JCOLUMN 
2 3 4 


IR0W=1 

IR0W=2 

IR0W=3 

IR0W=4 


***************** 
* * * * * 

* 1 * 5 * * * 

***************** 

# * # * * 

* 2 * 6 * * * 

* * **** ********** * 

* * * * * 

*3*7* * * 

***************** 

* * * * * 

*4*8* * * 

***************** 


IROW=l 

IR0W=2 

IR0W=3 

IR0W=4 


***************** 

* + * + * * * 

* 1 * 3 * 2 *' 4 * 

***************** 

* , * * * * * 

* 9 *11 *10 *12 * 

***************** 

* ♦ * * * * 

*7*5*8*6* 

***************** 

* + * * * * 

*15 *13 *16 *14 * 

***************** 


DISK FILES ON WHICH CHUNKS 

are stored 


ORDER IN WHICH CHUNKS ARE 
MULTIPLIED BY APPROPRIATE 
PORTION OF TRANSFORMED GREENS 
FUNCTION 


C0MM0N/ALLCOM/N,N02»N21»N04»N41»N34,NH,NHO2,NH21»NCH,NRHO,NHH, 

1 I2A,I2B,I3A»I3B,NH04 

COMMON/TRANCOM/RHOl (16,16*8) ,RH02 ( 16, 16,8) ,RHQ3 < 16, 16, B) , 

1 RH04(16, 16,8) ,HH(16, 16,9) 

COM-MON/HN21COM/HN21 (33,9) 
if C THE INITIAL CONDITIONS OVERLAY (PROGRAM INITIAL) OR PARTICLE ADVANCING 
C OVERLAY (PROGRAM STARS) STORES THE DENSITY CHUNKS OF IROW 1 AND 2 FOR 

I 


C JC0LUMN=1 ON DISK FILES 1 AND 2 RESPECTIVELY AND FOR JCOUiMN»2 ON DI&K FILES 
C 5 AND 6 RESPECTI VLEY. THE GETPHI OVERLAY REPLACES THE DENSITY ON THESE DISK 
r Fri fS utth the PDRRFSponOtNG VALUES OF POTENTIAL WHICH ARE THEN JSED IN THE 
C PARTICLE ADVANCTNG OVERLAY. THIS IS ACCOMPLISHED THROUGH CALLING SUBROUTINES 
C ANLX ( JCOLUMN) * ANLSYN < IROW > , AND SYNX ( JCOLUMN) AS DETAILED BELOW. 


C SUBROUTINE ANLX ( JCOLUMN) READ 5 RESPECTIVELY IROW 1 AND 2 FROM THE FOLLOWING 
r niSK FILES - 1 AND 2 FOR JCOLUMN=l AND 5 AND 6 FOR JC0LUMN=2 c IT THEN 
C PERFORMS A PERIODIC ANALYSIS IN THE X DIRECTION OVER JCOLUMN FOR I=1»N AND 
C WRITES THE RESULTS RESPECTIVELY FOR IROW 1*2*3* AND A ON THE FOLLOWING DISK 
C FILES - 1,2*3* AND 4 FOR JCOLUMN = l AND 5*6*7* AND 8 FOR JCOLUMN=H. 

CALL ANLX(l) 

C A L L ANLX ( 2 ) _ ^ . . r*r>i ■ ai i t mr 

C SUBROUTINE ANLSYN(IROW) READS RESPECTIVLEY JCOLUMN : 1 AND 2 FROM' THE FOLLOWING 
C DISK FILES - 1 AND S FOR IR0W=1,2 AND 6 FOR IROW=2,3 AND 7 FOR IR0W=3, AND 
C 4 AND 8 FOR IR0W=4. IT THEN PERFORMS A PERIODIC ANALYSIS IN THE Y DIRECTION 
C OVER IROW FOR J=1*N. FOR EACH CHUNK* IT THEN PERFORMS A PERIODIC ANALYSIS IN 
C THE Z DIRECTION FOR K=1 *NH« ELEMENT BY ELEMENT MULTIPLICATION WITH A 
C SIMILARLY SHAPED CHUNK OF THE TRANSFORMED GREENS FUNCTION AND THEN A PERIODIC 
C SYNTHESIS IN THE Z DIRECTION FOR K=1»NH. THE RESULT FOR K=l*NH/2 I5 T J H . E ^ C _ 

C PERIODICALLY synthesized IN the Y DIRECTION OVER IROW FOR J-l.N. THIS LAST 
C RESULT FOR JCOLUMN 1 AND 2 IS WRITTEN* RESPECTIVLEY ON THE FOLLOWING DISK 
r FTl FS - 1 AND 5 FOR IROW=l* 2 AND 6 FOR IROW^P* 3 AND 7 FOR IRQW=3* AND 4 AND 
C 8 FOR IROW=4. THE ORDER IN WHICH ANLSYN IS CALLED FOR IROW 1 THROUGH 4 
C MINIMIZES READING FROM DISK FILE 9 OF CHUNKS OF THE TRANSFORMED GREENS 

C FUNCTION AS MENTIONED ABOVE. 

CALL ANLSYN (1) 

CALL ANLSYN { 3 ) 

CALL ANLSYN (2) 

CALL ANLSYN(4) 

C SUBROUTINE SYNX (JCOLUMN) READS RESPECTIVELY IROW 1*2*3* AND 4 FROM T H E 
P Fm l OWING DI5K FILES - 1*2*3* AND 4 FOR JCOLUMN = l AND 5*6*7* AND 8 FOR 
C JCOLUMN=?, IT THEN PERFORM5 A PERIODIC SYNTHESIS IN THE X DIRECTION OVER 
C JCOLUMN FOR J=1*N. IT THEN WRITES THE RESULT RESPECTIVELY FOR IROW 1 AND 2 
C ON THE FOLLOWING DISK FILES - 1 AND 2 FOR JCOLUMN^l AND 5 AND 6 FOR JC0LUMN=2. 
CALL SYNX(l) 

CALL SYNX (2) 

return 

END 


SUBROUTINE ANLX(JCOLUMN) 
COMMON/ALLCOM/N*N02*N21»NOA*NA1* 


N3A*NH*NK02*NH21 *NCH*NRHO»NHH* 


1 

1 


I2A*I2B,I3A»I3B*NH0A 
COMMON/TRANCOM/RHOl (16*16*8) * 
RHOA (16*16*8) * HH (16*16*95 


RH02 (16*16*8) *RH03 (16*16*8) 


* 


COMMON Z(1025>, Y(1025) 
IF(JCOLUMN.EQ.2> GO TO 2 
READ ( 1 ) RHOl 


REWIND 1 
READ (2) RH02 


REWIND 2 
GO TO 3 

2 READ ( 5 ) RHOl 
REWIND 5 
READ (6) RH02 
REWIND 6 

3 CONTINUE 

CALL GETSET(3»I2A> 

DO 10 K=1*NH02 
DO 10 Js=l * NOA 

DO 5 1 = 1 * NOA 
Z ( I > =RH01 ( I * J»K) 

Z (NOA+I ) =RH02 ( I *J*K) 

Z (N02+ I ) =0 • 

5 Z (N3A+ I ) =0 • 

CALL FTRANS (3* I2A) 

DO 10 I=l*NOA 
RHOl (I»J»K)=Y(I> 

RH02 ( I * J*K ) =Y (NQA+ I ) 
RH03CI,J*K)=Y(N02+I) 

10 RHOA ( I * J *K ) =Y (N3A+ I ) 

IF(JC0LUMN.EQ.2> GO TO 12 

WRITE(l) RHOl 

REWIND 1 

WRITE (2) RH02 

REWIND 2 

WRITE (3) RH03 

REWIND 3 

WRITE (A) RHOA 

REWIND A 


\ 


RETURN . 
12 WRITE (5) 
REWIND S 
WRITE £ 6 > 
REWIND 6 
WRITE (7) 
REWIND 7 
WR ITE (8 ) 

REWIND 8 
RETURN 

END 


RH01 

RH02 

RH03 

RH04 


SUBROUTINE ANLSYN(IROW) 

C0MM0N/ALLC0M/N»NO2*N21 *N04*N41 *N34*NH*NH02*NH21 *NCH*NRHO*NHH* 


1 I2A» I2B* I3A* I3B*NH04 

COMMON/TRANCOM/RHOl (16*16*B) *RH02 ( 16* 16*8) *RH03 < 16* 16*8) * 
1 RH04 (16*16*8) *HH (16*16*9) 


COMMON/HN21COM/HN21 (33*9) 
COMMON Z ( 1025) * Y(1025) 


GO TO (l»2*3t4> IROW 


1 READ(l) RHOl 
REWIND 1 


READ (5) RH02 
REWIND 5 
GO TO 5 

2 READ(2) RHOl 
REWIND 2 
READ (6) RHQ2 


REWIND 6 


GO TO 5 

3 READ<3) RHOl 
REWIND 3 
READ (7) RH02 
REWIND 7 


GO TO 5 

4 READ (4) RHOl 
REWIND 4 

READ (8) RH02 
REWIND 8 

5 CALL GETSET(3*I2A) 

DO 10 K=1*NH02 

DO 10 I — 1 * NO 4- 
00 7 J=1 » N04 

Z ( J) “RHOl ( I * J*K ) 
Z(N04+J)=RH02(I*J*K) 
Z (N02+ J) =0 • 

7 Z (N34+J) =0 * 

CALL FTRANS ( 3 * I2A ) 

DO 10 J=1*N04 
RHOl (I* J*K)=Y(J) 

RH02 ( I » J»K) =Y (N04* J) 
RH03 ( I * J»K) -Y (N02+J) 


10 RH04(I,J»K)=Y(N34*J) 

GO TO (49,49*75,75) IROW 

49 READ (9) HH 

50 JCQLLJmN=1 

DO 70 1 = 1 * N04 
DO 70 J=1 » N04 
DO 52 K=1,NHD2 
Z<K)=RH01 (I« J*K> 

52 Z (NH02+K ) =0 • 

CAUL GETSET (3* I3A) 

CALL FTRANS (3* I3A) 
IF(IR0W.NE.3) GO TO 300 
IF(I.NE.l)GO TO 300 

LL=J 

GO TO 200 
54 DO 70 K=1,NH02 
70 RH01 (I,J»K)=Y(K> 

GO TO 100 

74 READ (9) HH 

75 JCOLUMN-2 

DO 95 1=1, N04 
DO 95 J=1,N04 
DO 77 K=1,NH02 
Z(K)=RH02(I, J,K) 

77 Z (NHD2+K ) =0 • 

CALL GET5ET(3,I3A) 

CALL FtRANS (3* I3A) 

IF (IROW. ME. 3) GO TO 300 

IF ( I »NE. 1 ) GO TO 300 

LL=N04+J 
GO TO 200 

79 DO 95 K=1,NH02 
95 RH02 (I*J»K)=Y (K) 

GO TO 125 

100 JC0LUMN=3 

DO 120 1=1, N04 
DO 120 J=1,N04 
DO 101 K=1*NH02 
Z(K)=RH03(I*J,K) 

101 Z(NH02+K)=0. 



CALL GETSET ( 3* I3A) 

CALL FTRANS { 3 • I3A) 

GO TO (103*105*107*115) IROW 
103 IFtJ.NE.DGQ TO 300 
LL=I 

GO TO 200 

105 IF ( J.NE.l) GO TO 300 
LL-N04+ I 

107 IF (I .N£. 1 • AND • J *NE • 1 ) GO TO 300 
IFd.EQ.l • AND. J.EQ • 1 ) GO TO 111 
IFtI.EQ.DGO TO 109 
LL=I 

GO TO 200 

109 LL=J 

GO TO 200 

111 LL=N21 

GO TO 200 

115 IF(J.NE.l) GO TO 300 
LL=N04+I 
GO TO 200 

117 DO 120 K=1*NH02 
120 RH03CI* J»K)=Y(K) 

| GO TO (74*74*400*390) IROW 

| 125 JC0LUMN=4 

I DO 145 1-1 *N04 

I DO 145 J=1*N04 

I DO 127 K=1*NH02 

Z (K) -RH04 ( I , J*K ) 

127 Z (NH02+K) s 0 • 

CALL GETSET C 3 * 1 3 A ) 

CALL FTRANS ( 3 * I3A) 

IF ( IR0W.NE.3) GO TO 300 
IF(I.NE.1)G0 TO 300 

LL=N04+ J 
GO TO 200 

129 DO 145 K=1*NH02 
145 RH04 ( I * J*K) =Y (K) 

GO TO (400*400*49*49) IROW 
200 DO 205 K=2,NH02 





Z(K)=Y(K)*HN21 (LL*K) 

205 Z (NH02+K ) S Y (NH02+K) tt HN2l (LL'*K> 
Z(1)=Y(1)*HN21 (LL ♦ 1 > 
Z(NH21)=Y(NH21)»HN21 (LL*NH2l> 

GO TO 310 

300 DO 305 K=2*NH02 

Z(K)=Y(K)*HH(I*J*K) 

305 Z (NH02 + K ) ~Y ( NH02+K > ^HH ( I * J*K) 
Z(1>=Y(1>*HH(I*J*1> 
Z(NH2l)=Y(NH21>*HHUtJ»NH2D 

310 CALL GETSET(4*I3A) 

CALL FTRANS(A»I3A) 

GO TO (54*79* 117* 129) JCOLUMN 
390 REWIND 9 
400 CALL GETSET(4*I2A) 

DO 410 K=:1*NH02 

DO 410 1-1 *N04 
DO 405 J=1*NQ4 
Z( J)=RH01 (I* J»K) 
Z(N04+J)=RH02(I»J»K) 
Z(N02+J)=RH03(I*J*K) 

405 Z (N34+J) =RH04 ( I * J*K) 

CALL FTRANS (4* I2A) 

DO 410 J-l *N04 
RHQ1 ( I * J*K) =Y ( J) 

410 RH02 (I*J»K)=Y (N04+J) 

GO TO (415*420*425*430) IROW 
415 WRITE ( 1 ) RHOl 
REWIND 1 
WRITE (5) RH02 
REWIND 5 

RETURN 

420 WRITE ( 2) RHOl 
REWIND 2 
WRITE (6) RH02 
REWIND 6 
RETURN 

425 WRITE (3) RHOl 

REWIND 3 
WRITE (7) RH02 



REWIND 7 
RETURN 

430 WRITE (4) RHOI 
REWIND 4 
WRITE (8) RH02 
REWIND 8 
RETURN 
END 


SUBROUTINE SYNX<JCOLUMN> 

COMMON/ALLCOM/N,NO2,N21.N04»N41,N34,NH»NH02,NH2i.NCH,NRH0,NHH, 
1 I2A.I2B.I3A, I3B.NH04 

COMMON/TRANCOM/RHOl <16,16,8>.RH02<16,16,B) ,RH03 ( 16* 16, 8 ) . 

1 RH04 (16.16.8) .HH (16.16, 9) 

COMMON Z(1025) *Y(1025) 

IF(JC0LUMN.EQ.2> GO TO 1 

READ ( I ) RHOl 

REMIND 1 

READ ( 2 ) RH02 

REMIND 2 

READ (3) RH03 

REWIND 3 

READ (4) RH04 

REWIND 4 
GO TO 4 

1 READ (5) RHOl 

REWIND 5 
READ (6) RH02 

REWIND 6 
READ(7> RH03 
REWIND 7 
READ (8) RH04 
REWIND 8 

4 CALL GEtSET (4* I2A) 

DO 10 K=1 .NH02 

DO 10 J=1,N04 
DO 5 1=1. N04 
Z(I)=RH01 (I, J.K) 

Z (N04+ I ) =RH02 ( I y J.K) 

Z (NQ2+I ) =RH03 ( I ♦ J.K) 

5 Z(N34+I)=RH04( I, J.K) 

CALL FTRANS ( 4. I2A) 

DO 10 1=1, N04 

RHOl (I*J,K)=Y(I) 

10 RH02(I,J,K)=Y(N04+I) 

IF(JC0LUMN.EQ.2) GO TO 12 

WRITEU) RHOl 
REWIND 1 
WRITE (2) RH02 


REWIND 2 
RETURN 
WRITE (5) 
REWIND S 
WRITE(6> 
REWIND 6 
RETURN 
END 


RHOl 

RH02 


o o 


■ OVERLAY (IFILE»3»0) 

C THlS P OVERLAY I SPECIFIES INITIAL c °SiDl TI °NS, GENERATES INITIAL ^TAR POSITIONS 
C AND VELOCITIES AND CREATES THt INITIAL DtNS ITY MESH. 

COMMON/ ALLC0M/N*N02*N21 *N04*N41 » N34,NH*NH02»NH21*NCH*NRHO»\JHH* 

* COMMON/ ADVC0M/N8R»NBS, Rl , XM»DT » DTE2 « N04M1 ? NOAM 2 » VXYMAX ♦ VZMAX I » 

1 CY*CYY»NPLOT*NPRINt* IN (2) ♦xmin*xmax*ymin»ymax»zmin, 

? ZMAX.RMIN.RMAX* VRMIN, VRMAX . VTMIN.VTMAX.VZMIN.VZMAX *XPP ' (3) * 

3 YP t ZP » YPP (3) .RPP ( 3) » VTP» VRP» VZP» IT APX. P I . MASKl * S? * JT * J s » NMK * 

4 NaS3*CXY»CZ»VZAVl *DR» ITEST . JTFILE* JSFILE*DDD*N04P1 
COMMON/ 1 NITCOM/MASSDd 6*4) 

DIMENSION PHI (32.32.8) 

DIMENSION XS(6)»RS{6) 

DIMENSION XPACK (20A8) .YPACK (2048) .ZPACK (2048) 

DIMENSION RPACK (2048) * TP ACK (204-8 ) 

DIMENSION E ( 16*4) * S IGMAVZ <16*4) *VTR < 16*4) .VRR (16*4) *TH( 16* 4) ♦ 
1 EQUIVALENCE 6 USU)?Xw"xlt2?!vX),(*S(3>,V),(XS(*),VY),<XS<5),ZI, 

1 EQUIVALENCE (RS(1) ,R) * (RS (2) *VR) * (RS (3) ♦THETA) * (RS(4) *VT>* 

1 (RS(5) .ZR) » (RS(6) »VZR) 

EQUIVALENCE (XPACK (I) .RPACK (1) ) * (YPACK(l) t TRACK (1) ) 

integer q.cy.cyy 

REAL MASSD 
CALL PSEUDO 

C IF THIS IS THE SECOND CALL TO THIS OVERLAY GO TO 70. 

IF ( ITEST.EQ.1 ) GO TO 70 
SET INITIAL VALUES 

SET TOTAL NUMBER OF STARS 
NBR=20480 

C SET NUMBER OF STARS PER READ OF STAR FILE 
NBS=2048 

c set initial radiusimust be ,le* n/4-3) 

RI“1 2 * 

c set initial standard deviation of z at zero radius 

c SET^ INITIAL MAXIMUM DEVIATION OF Z at ZERO RADIUSIMUST BE .LE. Nri/4-1) 

c SET X NUMBER OF TIME STEPS PER ROTATION 


' 

<1 


! DDD-50 

I c set scaled mass 

$ XM=3.55E6/NBR 

| C SET NUMBER OF POINTS PER PLOT 

| NP"NBR 

| C 5ET INITIAL MAXIMUM ABSOLUTE VALUE OF VX AND V Y 
■ | VXYMAX =3. , _ , 

4 c SET INITIAL MAXIMUM ABSOLUTE VALUE OF VZ 

I VZMAXI=.5 , _ 

| C SET TOTAL NUMBER OF TIME StEPS 

J. CY = 8 

I c SET NUMBER OF TIME STEPS BETWEEN PLOTS (MUST BE.GE.6) 

I NPL0T=6 

I C SET NUMBER OF TIME STEPS BETWEEN PRINTED DIAGNOSTICS 

I NPRINT = 2 _ £ . 

If C SET PLOTTING SPECIFICATIONS 

C SET PLOT TITLE 

IN ( 1 1 =1 OHJ • MILLER 
IN (21 =1 OH 3D- TEST 

C SET MAXIMUM AND MINIMUM VALUES TO BE PLOTTED 

XMIN=-8. 

XMAX=40. 

YMIN=-8. 

YMAX=40. 

ZMlN=-8. 

ZMAX=40. 

RMIN=0. 

RMAX-24 • 

VRMIN=-10. 

VRM AX= 1 0 • 

VTMIN=-10. 

VTMAX=10. 

VZMIN=-10. 

VZMAX=I0. 

C SET COORDINATE LABELS 

XPP ( 1 ) =1 OHX * 

XPP (2)=10HCYCLE= 

YP=I OH Y 
ZP=10H Z 
YPP ( I ) =1 OHY « 



YPP ( 2 ) =1 OHCYCLE= 

RPP ( 1 ) =1 OHR * 

RPP (2)=10HCYCLE= 

VTP = 1 OH VT 'ETA 

VRP“1 OH VR 
VZP=10H VZ 

C SET TAPE NUMBER FOR DDIPLT 

ITAPX=6LTAPE22 

C SET CONSTANTS 

PI=3. 1415926536 

MASK1 =077 777777770000000000 

S2=S0RT(2.) 

C INITIALIZE STAR TAPE NUMBERS 

JT=10 
JS = 1 1 
REWIND 10 
REWIND 11 
NMK=NP-NBS 
N8S3=3*NBS 


CXY=N04 

CZ=NH04+.5 

N04Ml=N04-l 

N04M2=N04-2 

DR=XM/2 

DR=DR. AND. MASK! 

RI2=RI*RI 

C END OF CONSTANTS 

C SET PHI TO ZERO FOR BUILD UP OF X-Y-Z MASS DENSITY 
DO 25 K=1,NH02 
DO 25 J=1,N02 
DO 25 1=1, N02 


25 PHIU,J,K)=0. 

C SET MASbD TO ZERO FOR BUILD UP OF AN AXlSYMMETRIC (R-Z) MASS DENSITY 
C WHICH WILL BE USED LATER IN INITIALIZING PARTICLE VELOCITIES. 

DO 30 IZ=1»NH04 


DO 30 IR=1,N04M1 
30 MASSD(IR,IZ)=0. 

C GENERATE RANDOM STAR POSITIONS WITH A RADIAL DENSITY DISTRIBUTION OF 
C SORT (1- (RADIUS/MAXRADIUS) *»2) AND A MAXWELLIAN_Z DENSITY DISTRIBUTION. BUILD 
C UP AN X-Y-Z MASS DENSITY IN PHI 


AND AN R-Z MASS DENSITY IN MASSD. TEMPORARILY 


o o 


C ASSIGN ZEROES AS THE STELLAR VELOCITY COMPONENTS. PACK THE POSITION AND 
C VELOCITY COMPONENTS INTO EACH WORD OF XPACK* YPAC<» AND ZPACK AND WRITE 
C THOSE ARRAYS ON DISK FILE 10. 

NS2=0 

X=URAN ( 72737 • ) 

35 DO 60 IS=1,NBS 
40 X=2.*(URAN(0.)-.S) 

Y=2.*(URAN<0.)-.5) 

Z7=2.*(URAN(0.)-.5) 

RR=X*X+Y*Y+Z7*Z7 

IF{RR.GT.l.) CO TO 40 

X=RI*X 

Y=RI*Y 

R2=X*X+Y*Y 

R=SQRt(R2) 

IF(IS.EQ.NBS) R={NS2+NBS)#RI/N8R 
IF(R.LT.l.E-5) X=0.0001 
THETA=ATAN2(-Y,X) 

X=R*COS (THETA) + CXY 

y=-R*S IN ( THET A ) +CX V 
ZR=Sqrt ( 1-R*R/RI2) 

ZZ=S2*ZR°SIGMAZ*SqRT C-ALoG ( 1 • -URAN ( 0 • ) ) ) 

Z=ZZ*COS (2.«-PI*URAN (0 . ’) / 

ZM=ZR*ZMAXI 

IF (ABS(Z) .GT.ZM) Z=SIGN(?M,Z) 

Z=Z+CZ 

vx=o. 

VY = 0. 

VZ=0. 

IR=R + 1 .5 

IZ=ABS (Z-CZ) +1.5 

BUILD UP R-Z DENSITY IN MASSD (DIVISION BY AREA OF MASS RINGS WILL BE 

PERFORMED LATER) . 

MASSD (IR»IZ)=MASSD(IR»IZ) +DR 
IX=X+ .5 
JY=Y+ .5 
KZ=Z+.5 

C BUILD UP X-Y-Z DENSITY IN LAST HALF OF EACH WORD OF PHI. 

CALL DENS(DP»PHI(IXf JY,KZ) ) 

CALL PACK (XS (1 ) ,XPACK CIS) ♦YPACKCIS) ,ZPACK(IS) ) 


60 CONTINUE 

WRITE < 10 > XPACK»YPACK*ZPACK 
NS2=NS2+NBS 

IF(NS2.LT.NBR) GO TO 35 
REWIND 10 

GO TO 800 _ 

c the second call to this overlay begins here. 

C READ THE POTENTIAL INTO THE pHI MESH FROM DISK FILES 1.2, 5» AND 6. 
70 READ(l) ( ( (PHI ( I, J,K) , I=1»N04) ,J=1,N04) »K=1,NH02) 

REWIND I 

READ (2) ( ( (PHI ( I , J»K ) , I=NO/+Pl »N02) » J=1 ,N04) , K=1,NH02> 

READ^S) 2 ( ( (PHI < I » J*K) ♦ 1=1 »N04) , J=N04P1 ,N02) »K = 1 »NH02) 

REWIND 5 

READ (6) ( ( (PHI ( If J»K) , I=N04P1 , N02) , J=N04P1 ,N02) *K=1*NHQ2> 

REWIND 6 

C SET LEAST SIGNIFICANT HALF OF EACH WORD OF PHI MESH EQUAL TO ZERO. 

DO 60 K=1,NH02 
DO 80 J=1 , N02 


DO 80 1 = 1 » N02 

80 PHI (I,J,K)=PHI (It J,K> . AND. MASK 1 
C AVERAGE THE Z COMPONENT OF GRAVITATIONAL 

c store it in the e array as a function of 

DO 130 1=1 ,N04M1 
DO 120 K=2,NH04 


FIELD ALONG FOUR RADIAL DIRECTIONS. 
R AND Z AND PRINT IT. 


KP=NH04+K 


KM=NH04+ 1-K 

EIPKP=2.* (PHI ( NO 4+ I » N04 , KP- 1 ) -PHI (N04+I, N04.KP) ) 
EIMKP=2 . * (PHI (N04-I ,NG4»KP-I)-PHI (N04-I * N04, KP) ) 

EJPKP=2.* (PHI (N04.N04+ I,KP-1)-PHI {N04,N04+I ,KP) ) 
EJMKP=2.* (PHI (N04»N04-I »KP-1 > -PHI ( N04 , N04- I ♦ KP ) ) 

EKP= ( E IPKP+EIMKP+EJPKP+EJMKP ) /4. 

EIPKM=2.* (PHI (N04+I »N04»KM+1 ) -PHI (N04+I ?N04»KM) ) 
EIMKM=2«* (PHI (N04-I ,N04»KM+1 ) -PHI (N04-I , N04t KM) J> 
EJPKM=2. fl - (PHI ( N04 « N04+ I »KM + 1 ) -PHI (N04, N04+I ,KM) ) 
EJMKM=2.* (PHI (N04.N04-I »KM+1 ) -PHI (N04» N04-I t KM) ) 
EKM=(EIPKM+EIMKM*EJPKM+EjMKM)/4. 

120 E (1+1 ,K)=(EKP+EKM)/2. 

130 E(I+1»1)=0. 

DO 140 K=2»NH04 






KP=NH04+K 

km=nho4+i-k 

(N04»N04*KP-1)-PHI (N04»N04,KP) ) 

EKM=2 # *> (PHI <N04,N04,KM+1 ) -PHI (N04, N04,KM) ) 

140 E ( 1 *K) = (EKP + EKM)/2* 

E (1,1)=0. 

PRIMT 910 

910 FORMAT ( 1H0 *Z COMPONENT OF FIELD*) 

PRINT 920 » { (E ( IR» IZ ) * I Z= 1 ♦NH04) * IR = 1 ,N04M1 ) 

920 FORMAT ( 1 H ,4E14.7) 

DIVIDE THE RING MASSES IN MASSq Br RING AREA TO STORE MASS DENSITY IN MASSD. 
DO 160 1 2 ~ 1 , NH04 
DO 150 IR=2,N04Ml 

150 MASSQ (IR,IZ)=MASSd(IR,IZ>/(4 6 *PIMIR-1) ) 

160 MASSD ( 1 * IZ) =MASSD ( 1 » IZ ) / < «■ 50*P1 ) 

COMPOTE A STANDARD DEVIATION OF THE Z VELOCITY WHICH WILL CREATE A PRESSURE 
GRADIENT IN Z JUSt SUFFICIENT TO BALANCE THE Z COMPONENT OF THE GRAVITATIONAL 
FIEL NH04Pi2nh04+1 STANDARD DEVIATION IN SIGMAVZ AS A FUNCTION OF R A ND Z. 

DO 180 IR = 1 , N04M1 
SIGMAVZ ( IR»NH04) =0 . 

DO 180 I IZ=2 , NH04 
IZ=NH04P1-I IZ 


180 


D0 G 190 Z IR=i!nO 4M1 GHAVZ(IR * IZ + 1> +MA5SD(IR,IZ + 1)ftE<IR ’ IZ + 1) 


DO 190 IZ=1,NH04 

IF (MASSD(IR*IZ> • EQ • 0 • ) Go TO 185 
IF (SlGMAVZ(IR»lZ) . LT • 0 . ) SIGMAVZ ( IR* IZ) =0. 

SIGMAyZ (IR* IZ) =SQRT (SIGMAV7 (IR» IZ) /MASSDl IR, IZ) ) 

GO TO 190 

185 SIGMAVZ ( IR* IZ) =0. 

190 CONTINUE 

?JSU L T C °r P0NE ^ T 0F THE GRAVITATIONAL FlELO ALONG FOUR RADIAL 

directions, store it in the e array as a function of r and z* 

DO 200 12=1 , NH04 
«P=NH04+IZ 
KM=NH04+1-IZ 
DO 200 IR=1 , N04M2 
IP=N04+ IR 
IM=N04-IR 


AND PRINT IT' 


EIPKP =PHI (IP-1 ,N04,KP) -PHI (IP+1,N04»KP) 

EIMKP sPHI(IM+l,N04»KP)-PHI(IM-l,N04»KP) 

EJPKP =PHI (N04. IP-1,KP)-PHT <N04, IP+1,KP) 

EJMKP =PHI <N04,IM+1,KP)-PHI (N04, IM-1 »KP ) 

EKP= (EIPKP+ EIMKP+ EJPKP+ EJMKPJ/4. 

EIPKM =PH I (IP-1,N04,KM)-PHI ( IP+1,N04,KM) 

EIMKM =PHI (IM+1,N04,KM)-PHI ( IM-1 ,N04»KM) 

EJPKM =PH I (N04* IP-1*KM)-PHI <!\I04* IP + 1 »KM) 

EJMKM =PHI (N04, IM+1 » KM ) -PH I < N04 , IM-1 ,KM) 

EKM= ( EIPKM+ EIMKM+ EJPKM+ EJMKM)/4. 

200 TH( IR+1, IZ)=(EKP+EKM)/2. 

DO 220 IR=2,N04Ml 
DO 210 IZ=2,NH04 

210 E(IR*IZ>-(TH(IR#IZ>+TH(IR*IZ-1) )/2. 

220 EUR* I)=3./2.*TH(IR»2)-.5*TH( IR,3) 

DO 230 IZ=1*NH04 
230 E (1»IZ)=0. 

PRINT 930 

930 FORMAT(1HO**R COMPONENT OF FIELDS) 

PRINT 920, ( (E( IR*IZ) * IZ=1»NH04) , IR=1*N04M1 ) 

C COMPUTE THE ANGULAR VELOCITY ANGV REQUIRED TO BALANCE THE RADIAL COMPONENT 
C OF THE GRAVITATIONAL FIELD AT HALF OF THE INITIAL RADIUS. 

IR02=R 1/2 • +1*5 

ANGV=SQRT( ABS ( E ( IR02 , 1 ) ) / ( IR02-1 ) ) 

C COMPUTE THE TIME OF ROTATION TROT » 

TROT=2.*PI/ANGV 

C USING the NUMBER OF TIME STEPS PER ROTATION DDD, COMPUTE THE TIME StEP DT. 
DT=TROT/DDD 
PRINT 242, XM 

242 FORMAT ( 1H0»3HXM~*E16*8) 

PRINT 244, DT 

244 FORMAT ( 1H ,3HDT=*E16.8) 

C TIME SCALE THE PARTICLE MASS DR AND SET lT$ LEAST SIGNIFICANT DIGITS EQUAL 
C TO ZERO. 

DT2=Dt**2 

DR=DP*DT2 

DR=DR.AND.MASK1 

C COMPUTE DtE2 WHICH IS LATER USED to APPLY A CENTRAL GRAVITATIONAL' FORCE TO 
C STARS WHICH ARE OUTSIDE OF THE MESh. 

DTE2=-XM*NBR*DT2 


c time scale the radial gravitational field e» the standard deviation of the 
c z component of velocity sigmavz and the r-z MASS DENSITY MAS5D» 

DO 246 I Z = 1 » NH04 
DO 246 IR=ltN04Ml 


E(IRtIZ)=E(IRtIZ)*DT2 

SIGMAVZ (IRtIZ)=SlGMAVZ(IR*lZ)*DT 

246 MASSDUR*IZ)=MASS0UR,IZ)*DT2 , rtC . lCr5AnT ,, 

C FOR EACH R-Z MESH POINT COMPUTE THE MINIMUM STANDARD DEVIATION OF THE RADIAL 
r ufi OCItY WHICH SATISFIES THE TOOMRE STABILITY CRITERION AND S T ORE IN SIGMAVR. 
C COMPUTE THE ANGULAR VELOCITY REQUIRED TO BALANCE THE RADIAL GRAVITATIONAL 
C FORCE AND STORE IN VTR. (COMMENTS WILL HENCEFORTH REFER TO THIS QUANTITY AS 

C THE ANGULAR VELOCITY OF THE COLD DISK.) 

DO 260 IZ=1»NH04 


DO 250 IR=2 t N04M2 

SIGMAVR <IR*IZ)=(E<IR+ltIZ>-E(IR-ltlZ>>/2.+E(IR»IZ>*3./(IR-l.) 

VTR (IRt IZ)=SQRT(E( IRt ID/ (IR-1 . ) > 

250 SIGMAVR < IR t IZ ) =3 . 36 *MASSq ( IR t IZ) *2.000/SQRT (ABS (SIGMAVR ( IRt IZ) ) ) 

SIGMAVR (1,IZ)=2.*SIGMAVR<2*IZ)-SIGMAVR<3,IZ> 

SIGMAVR (N04M It IZ)=SIGMAVR (N04M2’IZ) 

VTR (It IZ)=2.*VTR (2t IZ)-VTR (3t IZ) 

C FOR EAC^R-Z^ESH POIN^COMPUTE THE RATIO OF THE STANDARD DEVIATIONS OF THE 
C AZIMUTHAL AND RADIAL VELOCITIES AND STORE IN TH. 

DO 280 I Z= 1 1 NH04 
DO 270 IR=2 1 N04M2 

270 TH ( IRt IZ) =SQRT ( ABS ( 1 . + ( IR-1 . ) * (VTR ( IR+ 1* IZ) -VTR ( IR-1 t IZ ) > / 


1 ( 4 • *VTR ( IR t IZ ) ) ) ) 

TH(ltIZ)=l. 

280 TH(N04MltIZ)=TH(N04M2t IZ) 

C FOR EACH R-Z MESH POINT COMPUTE THE 
c square of the standard deviation of 


PRODUCT OF the MASS DENSITY AND THE 
THE RADIAL VELOCITY AND STORE IN VTR. 


c 

c 

c 

c 


DO 290 I Z = 1 1 NH04 
DO 290 IR=ltN04Ml 

290 VTRURtIZ)=MASSD(IRtIZ)*SIGMAVR(lRtIZ)**2 

FOR EACH R-Z MESH POINT COMPUTE THE DIFFERENCE IN THE SQUARES OF THE ANGULAR 
vFI OCItTES OF THE COLD AND WARM DISKS AND STORE IN VRR. THIS WILL NOt INCLUDE 
A TERM INVOLVING THE AVERAGE OF THE PRODUCT OF THE AXIAL AND RADIAL VELOCITIE 

WHICH WILL BE COMPUTED LATER. 


DO 310 IZ = 1 1 NH04 
DO 300 IR=2 1 N04M2 


IF(MASSD(IR,IZ).EQ.O.) GO TO 295 

VRR<IR,IZ)=(VTR< IR*1*IZ>-VTR(IR-1*IZ) >/( ( IR-1 . ) »MASSq < IR* IZ > *2. > + 

1 VTR<IR,IZ)*<l.-TH(IR,lZ)**2)/< < IR-1 . ) * ( IR-1 . ) *MASSD ( IR, IZ) ) 

GO TO 300 

295 VRR < IR, IZ) =0. 

300 CONTINUE 

VRR (1 . IZ)=2.*VRR (2* IZJ-VRR (3* IZ) 

IF (MASSD(3* IZ) .EQ.O.) VRR ( 1 * I Z ) =VRR (2, IZ) 

IF (MAS$D ( 1 * IZ) .EQ.O.) VRR ( 1 » IZ ) =0 . 

VRR (N04M1 , IZ)=VRR (N04M2, IZ) 

IF (MASSO (N04M1 * IZ) .EQ.O.) VRR (NOAM 1 * IZ) =0 . 

310 CONTINUE 

SET TO ZERO THE R-Z MESH VRVZ IN WHICH WILL BE SUMMED THE PRODUCT OF THE 
AXIAL AND RADIAL VELOCITIES. 

DO 320 IZ=1,NH04 
DO 320 IR=1 *NOAMl 

320 VRVZ ( IR » IZ ) =0 , 

READ THE PACKED SyAR POSITION AND VELOCITY COMPONENTS FROM TAPE 10. UNPACK 
THE POSITION AND VELOCITY COMPONENTS. DO A BILINEAR INTERPOLATION OVER THE 
R-Z MESH POINTS OF THE STANDARD DEV 7 AT ION OF THE RADIAL VELOCITY IN SlGMAyR 

and store the value in vR. do a bilinear interpolation of the ratio of the 

standard deviations of the azimuthal and radial velocities contained in th 

and store the value in vt. do a bilinear interpolation of the standard 

deviation of the axial velocity and store in vz. using these interpolated 
values of the Standard deviations of the velocities in the radial azimuthal 

AND AXIAL DIRECTIONS, GIVE THE PARTICLES MAXWELLIAN VELOCITY DISTRIBUTIONS 
IN THESE DIRECTIONS AND StORF IN vR* VT, ANO VZ. IN THE VRVZ MESh SUM THE 
PRODUCT OF THE RADIAL AND AXIAL VELOCITIES. COMPUTE VX AND VY FROM vR AND VT- 
PACK X AND VX, Y AND VY, AND Z AND VZ AND STORE ON DISK FILE 11. 

X=URAN ( 123A567 . ) 

VZAV 1=0. 

N S 2 = 0 

330 READ (10) XPACK,YPACK»ZPACK 
DO 360 IS=1,NBS 

CALL UNPACK (XS(1) ,XPACK ( IS) , YPACKT IS) ,ZPACK (IS) } 

XC='X~CXY 

YC=Y-CXY 

R2=XC*XC+YC*YC 

R=SQRt(R2) 

IR=R 


DRR=R-IR 
I R= I R + 1 
ZC=A85 (Z-CZ) 

IZ=ZC 

DZZ=ZC-IZ 
IZ= I Z + 1 

D11=1-DZZ-DRR+DRR*DZZ 

D12=DZZ*tl-DRR) 

D21=DRR* (1-DZZ) 

D22=DRR*DZZ 

IF{MASSO(IR»IZ) .EQ.O.) 011=0. 

IF (MASSD(IR*IZ+1) .EQ.O.) 012=0. 

IF(MASSD(IR+l»IZ).EQ.O.) 021=0. 

IF (MASSD(IR+1*IZ + 1) • tQ . 0 . ) D22=0. 

DSUM=Q11+D12+D21+D22 

VR=Dll«SIGMAVR(IR»IZ)+0l2 i} SIG^AvR(IR»IZ + l)+D21*5IGMAvR(IR + l*I2) + 

1 D22^5IGMAVR(IR+1»IZ+1) 

VR=VR/DSUm 

VT=011»TH(IR« IZ) +D12*TH<IR,IZ+1> +D21*TH { IR+1 * IZ) +U22*TH (I R ♦ 1 , IZ + 1) 
VT=VT/DSUM 

VZ=Dll*SlGMAVZ<IR,IZ)*0l2*SIiiMAVZ(lR*IZ*l)*D21*SIGMAVZ(IR*l*IZ>* 

1 022#SIGMAVZ(IR+1*IZ+1) 

VZ=VZ/DSUM 

VR=VR*SQRT <-ALOG< 1 • —UR AN { 0 . ) ) ) *^2 
THET1=2.*PI*URAN(0.) 

VT=VT*VR*SlN<THETl ) 

VR=VR»COS (THET1 ) 

THETA=ATAN2 (-YC.XC) 

VX=VR-«'COS (THETA) ~VT«- S I iN (tHETA) 

IFfABStvX) .GT.VXYMAX) VX=S IGN ( VXYMAX * VX ) 

V Y=-VR»S IN (THETA )-VT*COS (THETA) 

IF (ABS { VY) .GT.VXYMAX) VY=S I GN ( VXYMAX » VY > 
VZ=VZ*SQRT(-ALOG(1.“URaN(0.) ) )*52 
VZ=VZ»COS(2.*PI»URANCO.) ) 

IF ( ABS [ VZ ) .GT.VZMAXI) VZ=SJGN CVZMAXI»VZ) 

VZAVl=VZAVl+VZ 
IZ=ZC+1 .5 
IR=R+1 .5 

VRVZ (IR»IZ)=VRVZ fIR*IZ) +VR^VZ 
360 CALL PACK(XSd) t XPACK{lS) ,YPACK(IS) ,ZPACK(I5) ) 


WRITE (11) XPACK » YPACK » 2PACK 
NS2=NS2+NBS 

IF(NS2.LT.NBR) 60 TO 330 
REWIND 10 

REWIND 11 

VZAV1=VZAV1/N8R 

C ADD THE AXIAL PARTIAL DERIVATIVE OF THE AVERAGE VALUE OF THE PRODUCT OF THE 
C RADIAL AND AXIAL VELOCITIES TO WHAT IS ALREADY IN THE VRR MESH TO FORM THE 
C DIFFERENCE IN THE SQUARES OF THE BALANCED ANGULAR VELOCITIES OF THE COLD 
C AND WARM DISKS. 

NH0AM1=NH04-1 
DO 365 IR=2 ♦ NOAM1 
DO 3B3 IZ=2,NH0AM1 
IF(MASSD(IR.IZ> .EQ.O.) GO TO 383 

VRR ( iRt IZ)=VRR ( IR, IZ ) + < VRVZ < I R * IZ+ 1 ) -VRVZ ( IR ♦ I Z-l ) )*DR/ 

1 (2.«MIR-1)*MaS5D(IR*IZ) ) 

383 CONTINUE 

VRR(lR»l)=2e*VRR< IR,2)~VRR (IR»3) 

IF (MASSD ( IR* 3) • ECt • 0 • > VRR ( IR » 1 ) =VRR ( IR *2 ) 

IF (MASSD ( IR« I ) .EQ.O.) VRR(IR»1)=0. 

VRR (IR.NHOA) =VRR (IR»NH04M1 ) 

IF (MASSD (IR.NHO^ .EQ.O.) VRR { IR ♦ NHOA ) =0 . 

385 CONTINUE 

VRR(lfl)=2.*VRR(2$I)-VRR<3»l) 

IF (MASSD (3, 1 ) .EQ.O.) VRR ( 1 * 1 ) =VRR <2t 1 > 

IF(MASSD(1*1) .EQ.O.) VRR(1,1)=0. 

C ENCODE CYCLE NUMBER CYY In PREPARATION FOR PLOTTING. 

ENCODE (10,388,XPP(3) ) CYY 
ENCODE (10.388.YPP(3) ) CYY 
ENCODE (I CW388»RPP(3) ) CYY 
388 FORMAt(IIO) 

C READ THE STAR POSITIONS AND VELOCITIES from DISK FILE 11 AND UNPACK. DO A 
C BILINEAR INTERPOLATION OVER R-7 OF THE RADIAL GRAVITATIONAL FIELD IN E AND 

c store in er. divide by the radius to obtain the square of the balanced 
c angular velocity of the cold disk and store in angv. do a bilinear 
C INTERPOLATION OF VRR AND SUBTRACT FROM ANGV and TAKE Sqrt to yield average 
C BALANCED ANGULAR VELOCITY AT THAT POINT. ADD THE RESOLVED COMPONENT OF THIS 
C VALUE TO VX AND VY. SUBTRACT THE AVERAGE VAL l 'E OF VZ FROM V Z. TIME CENTER 
C THE POSITIONS. COMPUTE A NEW DENSITY IN THE l ,ST HALF OF EACH WORD OF THE PHI 
C MESH. PACK THE STAR POSITIONS AND VELOCITIES AND WRITE ON DISK FILE 10. MAKE 
C AN X-Y PLOT OF STAR POSITIONS. 


NS2 = 0 

390 READ (11) XPACK»YPACK»ZPACK 
DO 420 IS=1,NBS 

CALL UNPACK (XS{1) ,XPACK(IS) ,YPACK(IS) ,2PACK{IS) ) 

XC=X-CXY 

YC=T-CXY 

R=50Rt (XC*XC+YC*YC> 

IR-R 

DRR=R-IR 

IR=Ik+1 

ZC=AB5 (Z-CZ) 

IZ=ZC 

DZZ=ZC-IZ 

IZ=IZ+1 

D11-1-DZZ-DRR+DRR*DZZ 

D12=DZZ*(1-DRR) 

D21=DRR»(1-DZZ> 

D22=DRR*DZZ 

ER= D1 1*E ( IR» IZ) +D12*E < IR* IZ+1 ) +D21*E ( IR + 1 » IZ ) +D22*E < IR + 1 » IZ + 1 ) 
ER=ABS { ER ) 

RR = R 

IFtRR.LT. 0.11) RR=RR+0.05 
ANGV=ER/RR 

IF{MASSCMIR*IZ) .EQ.O.) Dll=0. 

IF(MASSD(IRf IZ+1) .EQ.O.) Dl2=0. 

IF (NASSD ( IR+1 » IZ) .EQ.O. ) D21=0. 

IF (MASSD(IR+1»IZ+1) .EQ.O.) D22=0. 

DSUM=D11+D12+D21+D22 

VP=D11*VRR <IR» IZ) +D12*VRR(IRtIZ+l) +D21*VRR ( IR+1 * IZ) + 

1 D22*VRR(1R+1»IZ+1) 

VP=VP/DbUM 
ANGV=AMGV-ABS (VP) 

IFtANGV.LT.O.) ANGV=0. 

ANGV=SQRT(ANGV) 

vy=vy-angv*xc 

VX=VX+ANGV°YC 

VZ=VZ-VZAVI 

THETA=ATAN2{-YC»XC) 

X=R*C0S{ANGV/2.+THETA) + CXY 


Y=-R*SIN (ANGV/2.+THETA) + CXY 

lX = X+.5 

JY=Y+.5 

KZ=Z+.5 

CALL DEM5(DR*PHI ( IX» JY,KZ) ) 

420 CALL PACK (XS ( 1 ) ,XPACK ( IS) , YPACK ( IS) ,ZPACK ( IS) ) 

Q = 0 

IF(NSa.EQ.NMK) Q = 1 

CALL ddiplt(Q,in*nbs,xpack,ypack,xmin,xmax,ymin,ymax, 

1 3*XPP*l*YP,13f ITAPX) 

WRITEQO) XPACK»YPACK*ZPACK 
NJS2=NS2+N8S 

IF(NS2.LT.NBR) GO TO 390 
REMIND 10 

C READ R THE N P0SI T I0NS and VELOCITIES FROM DISK FILE 10 AND MAKE AN X-Z PLOT OF 
c star positions. 

NS2 = 0 

430 READ ( 1 0 ) XPACKtYPACK.ZPACK 


Q=Q 

IFtNS2.EQ.NMK) Q = 1 

CALL DDIPLT (Q » IN » NBS ♦ XPACK , ZP ACK t XMIN » XMAX » ZM IN » ZMAX* 
1 3»XPP* 1»ZP»13»ITAPX) 

NS2=NS2+NBS 

IF (NS2.LT. NBR) GO TO 430 
REWIND 10 

READ FROM DISK FILE 10 AND MAKE A Y-Z PLOT 
X AND Y POSITIONS AND VELOCITIES TO RADIAL 
VELOCITIES AND WRITE THEM ON DISK FILE 12. 

NS2 = 0 

450 RE AD ( 1 0 ) XPACK »YPACK* ZP ACK 


of star positions, convert the 

AND AZIMUTHAL POSITIONS AND 


Q=0 

CALL^DDIPLT^Q. IN,NBS,YPACK,ZPACK,YMIN*YMAX,ZMIN,ZMAX,3,YPP* 

1 1,ZP»13, ITAPX) 

DO 480 IS=1,NBS 

CALL UNPACK (XS( 1) , XPACK ( I S ) ,YPACK (IS) ,ZPACK ( IS) ) 

XC=X-CXY 

YC-Y-CXY 

R=SQRt (XC»XC+YC*YC) 


TH£TA=ATAN2(-YC*XC) 

VR= (XC^VX+YC^VY) /R 

VT=-(XC«VY-YC*VX) /R 
ZR=Z 

VZR=VZ 

480 CALL PACK(RSd) ,RPACK(lS) ,TPACK<IS) ,ZPACK(IS) ) 

WRITE { 12) RPACK»TPACK, ZPACK 
NS2=.MS2 + NBS 

IF(N52.LT.NBR) GO TO 450 
REWIND 10 
REWIND 12 

C READ FROM DISK FILE 12 AND MAKE A PLOT OF RADIAL VELOCITY VS RADIUS. 
N|S2=0 

490 READ (12) RPACK , TPACK » ZPACK 
DO 510 15=1 *NBS 
510 TPACK ( IS)=RPACK (IS) 

CALL SHIFT2(RPACK(1) ,NBS) 

0 = 0 

IF (NS2.EQ.NMK) 0=1 

CALL DDIPLT (Q » IN,NBS»TPACK,RPACK»RMIN,RMAX» vRMINt VRMAX,3»RPP» 

1 1 , VRP ,13, ITAPX ) 

NS2=NS2+NBS 

IF (NS2.LT »NBR) GO TO 490 
REWIND 12 

C READ FROM DISK FILE 12 AND MAKE A PLOT OF AZIMUTHAL VELOCITY VS RADIUS. 
NS2=0 

520 READ (12) RPACK, TPACK»ZPACK 
CALL SHI FT 2 (TPACK ( 1 ) tNBS) 

Q=0 

CALLED IPLdo! IN,NBS,RPACK,TPACK,RMIN,RMAX,VTMIN,VTMAX,3f RPP, 

1 1 , VTP, l3» ITAPX) 

NS2=NS2+NBS 

IF(NS2.LT.NBR) GO TO 520 
' REWIND 12 

C READ FROM DISK FILE 12 AND MAKE A PLOT OF AXIAL VELOCITY V S RADIUS. 
NS2=0 

540 READ (12) RPACK, TPACK, ZPACK 
CALL SHIFT2(ZPACK (1) *NBS) 

Q=0 


IF (Ni2.EQ.NMK ) Q=1 

CALL DO IPLT (Q , IN, NBS ,RPACK , ZPACK , RM IN,RM Ax , VZM IN , VZMAX , 3 »RPP » 

1 1 »VZP» 13* ITAPX) 

NS2=NS2+NBS 

IF(NS2.LT.NBR) GO TO 5A0 
REWIND 12 

C SHIFT THE DENSITY TO THE MOSt SIGNIFICANT HALF OF EACH WORD OF THE pHI MESH. 
800 CALL SHIFT2(PHI (I,l» 1) *NCH) 

C WRITE the DENSITY FROM THE PHI MESH ONTO DISK FILES 1,2,5, AND 6. 

WRITE(l) ( ( (PHI (I, J,K) , 1=1, N04) t J=1,N04) ,K=1,NH02) 

REWIND 1 

WRITE (2) U (PHI ( I * J*K ) , I=N04P1,N02) , J=1»N04) »K=1»NH02> 

REWIND 2 

WRITE (5) ( ( (PHi ( I 9 J»K) *lssl,N04) , J=N04P1 ,N02) *K=1,NH02) 

REWIND 5 

WRITE (6) ( ( (PHI ( I » J,K ) * I=N04PI,N02) , J=N04P1 , N02) ,K = 1 ,NH02) 

REWIND 6 

RETURN 

END 



0VERLAY(SFILE,4*0) 

program stars 

C THIS OVERLAY USES THE GRAVITATIONAL POTENTIAL' COMPUTED IN THE GETPHI OVERLAY 
C TO COMPUTE NEW STAR POSITIONS AND VELOCITIES AND TO COMPUTE A DENSITY MESH 
C USED IN THE GETPHI OVERLAY* IT ALSO PRINTS OUT CERTAIN DIAGNOSTICS EVERY 
C CYCLE AND PERIODICALLY PRINTS OUT MORE EXTENSIVE DIAGNOSTICS AND MAKES PLOTS, 
COMM ON/ ALL CO M/N *N02*N21*N04*N41»N34*NH»NH02»NH21 *NCH*NRHO*NHH* 

I I2A*I2B*I3A,I3B«NH04 

COMMON/ A DVCOM /N8R »NBS ,RI ,XM*DT*DTE2*N04Ml*N04M2*VXYMAX*VZMAXl« 

1 CY*CYY,NPL0T*NPRINT*IN(2)*XMIN,XMAX*YMIN,YMAX,ZMIN, 

2 ZMAX»RMIN*RMAX*VRMIN*VRMAX*VTMIN*VTMAX*VZMIN*VZMAX»XPP (3) * 

3 YP*ZP»YPP(3) *RPP (3) ,VTP.VRP.VZP.ITAPX,PI,MASK1,S2,JT * JS* NMK* 

4 NBS3*CXY*CZ*VZAVl *DR» I TEST* JTFILE* JS FILE* ODD* N04P1 
DIMENSION PHI (32*32*8) ,MASSDR(16) *VR2AVGR(I6) »VRAVGR(I6) * 

1 VT2AVGR (16) *VTAVGR ( 16) *VZ2AVGR (16) *VZAVGR (16) »MASSDZ(8) * 

2 VZ2AVGZ(B) *VZAVGZ(B) 

DIMENSION XS(6) 

DIMENSION XPACK(2048)*YPACK(2048) *ZPACK(2048) 

EQUIVALENCE (XS(1) ,X) » (XS(2) *VX) *(XS(3) *Y) * (XS(4) *VY) * (XS(5) *Z) * 

1 (XS(6)*VZ> 

integer Q,CY*CYY*DDM 
REAL KE*MASSDR,MASSDZ 


f 

A 

I 
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•A 


CALL PSEUDO 
PRINT 20 

20 FORMAT (40H STAR NUMBER X VX » 

1 56H Y VY Z VZ 

2 R / ) 

C READ THE POTENTIAL INTO THE PHI MESH FROM DISK FILES 1*2*5* AND 6. 
READ ( 1 ) ( ( (PHI (I»J*K) *I=1*N04) *J=1*N04) *K=1,NH02) 

REWIND 1 

READ (2) ( ( (PHI ( I * J*K) * I=N04P1 »N02) * J=1*N04) ,Ki=1*NH02) 

REWIND 2 

READ (5) ( ( (PHI (I* J*K) *I=l*N04> *J=N04P1»N02) ,K=1*NH02) 

REWIND 5 

READ (6) ( ( (PHI ( I » J*K) * I=N04P1 »N02) * J-N04P1 *N02) *K-1 *NH02) 

RFWIND 6 

C SET TO ZERO THE LEAST SIGNIFICANT HALF OF EACH WORD OF THE PHI MESH 


. 1; 

DO 

60 

K s 

1 

* NH02 


DO 

60 

J= 

1 

♦ N02 

mm 

DO 

60 

1 = 

1 

♦ N02 

H 






m 1 







» 




f 

f 60 PHI (I,J,K)=PHI <I»J»K> .AND.MASK1 

|C DETERMINE WHICH PLOTS AND/OR PRINTED DIAGNOSTICS (IF ANY) ARE TO BE MADE THIS 
C CYCLE. 

t iprint=i 

IF(CYY-CYY/NPRlNT*NPRINT.EQ.O. .OR.CYY.EQ.l) IPRINT=2 
I IPL0T=1 

4 IF<CYY+5~(CYY+S)/NPLOT*NPLOT.EQ.O.) IPL0T=2 

j IF(CYY + 4“(CYY + 4)/NPLOT«-NPLOT.EQ.O.) IPL0T=3 

it IF(CYY+3-(CYY+3)/NPLOT*NPLOT.EQ.O.) IPL0T=4 

I IF(CYY+2-(CYY+2)/NPLOT*NPLOT.EQ.O.) IPL0T=5 

'] IF(CYY + 1"(CYY*1) /NPLOT*NPLOT .EQ.O.) I PLOT =6 

if IF(CYY-CYY/NPLOT*NPLOT.EQ.O.) IPLOT=7 

c IF NECESSARY ENCODE THE CYCLE NUMBER CYY INTO PLOT AXES LABELS. 

IF(IPL0T*EQ.2.0R • I PLOT • EQ • 3 ) ENCODE ( 1 0 *65 *XPP( 3 ) ) CYY 
I 65 FORMAT (I 10) 

:j IFUPLOT.EQ.A) ENCODE(10»65*YPP(3) ) CYY 

| IFUPL0T.GE.5) ENC0DE(l0 f 65*RPP(3) > CYY 

:|C IF LONG DIAGNOSTIC PRINTING iS TO BE DONE THIS CYCLE SET CONSTANT R MESHES 
C AND CONSTANT Z MESHES TO ZERO. 

IF{ IPRINT.EQ. I > GO TO 90 

1 DO 70 IR=l»NOA 

f MASSDR(IR).“0. 

VR2AVGR ( IR) =0. 

| VRAVGR ( IR ) =0 • 

VT2AVGR ( IR ) =0 • 

VTAVGR(IR)=0. 


VZ2AVGR ( IR ) =0 • 
70 VZAVGR(IR)=0. 

DO 80 KZ=1*NH02 


!c 

C 

3 c 

r|C 

ic 

c 


MASSDZ(KZ)=0. 

VZ2AVGZ (KZ)=0. 

80 VZAVGZ(KZ)-0. 

READ THE PACKED STAR POSITIONS AND VELOCITIES FROM DISk FILE UT • IE LONG 
DIAGNOSTIC PRINTING IS NOT TO BE DONE THIS CYCLE* THE STARS ARE PROCESSED 
THROUGH A FIRST DO LOOP ENDING WITH LINE 150. IF LONG DIAGNOSTIC PRINTING 
IS TO BE DONE THIS CYCLE THE SjARS ARE PROCESSED THROUGH A SECOND 00 LOOP 
ENDING WITH LINE 210. BOTH LOOPS TAKE THE GRADIENT OF THE POTENTIAL IN THE 
PHI MESH AND USE IT TO COMPUTE NEW STAR POSITIONS AND VELOCITIES. ADVANCE 
THE PARTICLES WHICH ARE OUtSiDE OF A CYLINDER TANGENT TO THE PHI MESH BY 
ASSUMING THAT THE TOTAL MASS OF THE GALAXY RESIDES AT ITS CENTER. SUM THE 


SO 

Ik 


f'-SL 






c number of particles outside of this cylinder, compute the DENSITY and store 
c it in the least significant half of each word of the phi mesh, pack the 
C NEW STAR POSITIONS and velocities and write them ON DISK file js. DEPENDING 
c ON THE VALUE of IPLOT (WHICH DEPENDS ON CYCLE NUMBER AND PLOTTING FREQUENCY) 
C ONE OF Six PLOTS MAY BE MADE. 


KE=0. 
PEIN=0 * 
PEOUT=Q. 

90 NUMBER=1 
DDM=0 


NS2=0 

95 READ(JT) xpack»ypack*zpack 
IFUPRINT.EQ.2) GO TO 155 

BEGINNING OF 5jAR ADVANCING LOOP FOR NO LONG DIAGNOSjIC PRINTING 
CALL 5 UNPACK (XS a ) ,XPACK < IS ) ,YPAC<( IS > *ZPACK < IS) ) 


XC=X-N04 

YC=Y-N04 

R=:SQRj(XC*XC+YC*YC) 

I X=X+ • 5 
JY = Y+ .5 
KZ=Z+.5 
JR=R + X • 5 

IF (KZ.LE. 1 ) GO TO 110 
IF (KZ.GE.NH02 ) GO TO 110 
IF (IR.GE.N04) GO TO 110 
VX=VX+PHKIX*1 tJY,KZ)-PHI (IX-lt JY,KZ> 
VY=VY + PHI(IX,JY*1,KZ)-PHI (IX,JY-1,KZ) 
VZ=VZ*(Z-KZ+.5)' n ’(PHI<IX*jY«KZ + l) -PHI ( IX* JY*KZ) ) 
1 *(KZ~Z+.5)MPHI(IX*JY»KZ)~PHI(IX*JY» KZ~1 ) ) 

GO TO 120 
110 ZC=Z-CZ 

R2=XC*XC+YC*YC+ZC*ZC 

R32=R2*SqRt(R2) 

E=DTE2/R32 

VX=VX*E*XC 

VY=VY+E*YC 


ooonoonn onoooon 


VZ=VZ+E*ZC 
120 X=X+VX 
Y=Y*VY 
Z=Z*V Z 
XC=X-N04 


140 

150 


YC=Y“N04 

R=SQRt(XC*aC+YC*YC) 

IX=X*.5 

UY=Y* ,5 

KZ=Z+.5 


IR=R+ 1 .5 

IF(KZ.LT.l) GO TO 140 
IFIKZ.6T.NH02) 60 TO 140 
IF (IR.6T.N04) 60 TO 140 
CALL DEN5(DR*PHI(IXtJY,KZ>) 
60 TO 150 


DDM=DDM+1 

CALL PACKtXSfl) ,XPACK (IS) ,YPACK(IS) ,ZPACK.(IS) ) 

GO TO 220 

BEGINNING OF STAR ADVANCING LOOP USED FOR- LONG DIAGNOSTIC PRINTING 

FOLLOWING MESHES SUM THE DESIGNATED QUANTITY FOR A PARTICULAR MASS 
VR2AVGR(IR> - RADIAL VELOCITY SQUARED 
VRAVGR(TR) - RADIAL VELOCITY 
VT2AVGR { IR ) - AZlMUTHAU VELOCITY SQUARED 
VTAVGR(IR) - AZIMUTHAL VELOCITY 
VZ2AVGR ( IR ) - AXIAL VELOCITY SQUARED 
VZAVGR(IR) - AXIAL VELOCITY 
MASSDR(IR) - NUMBER OF PARTICLES 

the following meshes sum the designated quantity for a particular 


MASS LAYER - 

MASSDZiKZ) - number of particles 
VZ2AVG2 ( KZ ) - AXIAL VELOCITY SQUARED 


VZAVGZ(KZ) - AXIAL VELOCITY 
THE LOOP ALSO SUMS THE POTENTIAL 


AND KINETIC ENERGIES 


DO 210 iSsltNBS 

CALL UNPACK(XS(1) ,XPACK(IS) ,YRACK<IS) »ZPACK<IS> > 
XC=X-N04 


YC=Y-N04 


R=SQRT(XC*XC+YC»YC> 
IX=X* .5 


0 

l 


JY=Y + 
KZ=Z + 


.5 

.5 


. THE 
RING - 


AXIAL 


155 


IR=R+ 1 ,5 

IF(KZ.LE.l) GO TO 160 

IF (KZ.GE.NH02) Gu TO 160 
IF< IR.GE.N04) GO TO 160 

VX=VX+PHI {IX+1,JY,KZ)-PHI UX-1, JY,KZ) 

V/Y=VY + PHI (IX»JY + 1,KZ)-PHI (IX»JY-1,KZ) 

VZ=VZ+ (Z-KZ*.5)*(PHI ( IX*JY,KZ+1)-PHI (IX, JY,KZ> ) 
1 ♦ (KZ-Z*.5)*(PHI(IX,JY,KZ)-PHI (IX*JY,KZ-1> ) 

GO TO 170 
160 ZC=Z-CZ 

R2=XC*XC+YC*YC+ZC*ZC 

R32=R2?SQRT(R2) 

E=DTE2/R32 

VX=VX+E*XC 

vy=vy+e*yc 

VZ=VZ+E»ZC 
170 X=X+VX 
Y=Y ♦ VY 
Z-Z*VZ 
XC=X-N04 
YC=Y-N04 

R=SQRt (XC*XC+YC*YC) 

IX=X+ . 5 
JY=Y+ , 5 
KZ=Z + *5 

IR-R+1 .5 

IF (KZ.LT • 1 ) GO TO 200 
IF (KZ *GT #NH02 ) GO TO 200 

IF ( I R i* GT «N04 ) GO TO 200 
CALL DENS(DR»PHI (IX,JY*KZ) ) 
MASSDR(IR)=MASSDR(IR)*1. 

IF (R.EQ.O. ) GO TO 180 
VR=(XC»VX*YC*VY)/R 
VT=(XC*VY-YC*VX)/R 
GO TO 190 

180 VR=SQRT(VX»VX+VY*VY) 

VT=0. 

190 VR2AVGR ( IR) =VR2AVGR ( IR) +VR*VR 
VRAVGR ( IR) =VRAVGR (IR) +VR 


n- n t ran 1 II l i t- i t . r i-in iM j i i >n— rm ii-rf - 
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VT2AVGR ( IR ) =VT2AVGR ( IR ) +VT»VT 
VTAVGR ( IR ) =VT AVGR ( IR> +VT 
VZ2=VZ*VZ 

VZ2AVGR<IR)=VZ2AVGR(IR>+VZ2 
VZAVGR { IR) =VZAVGR ( IR) + VZ 
MASSDZ(KZ)=MASSDZ{KZ) +1. 

VZ2AVGZ(KZ)=VZ2AVGZ(KZ) *VZ2 
VZAVGZ (KZ)=VZAVGZ (KZ) +VZ 
PEIN=PEIN*PHl < IX* JY*KZ) 

GO TO 205 
200 DDM=DOM+l 
ZC=Z-CZ 

RR=SQRT <XC*XC+YC*YC+ZC*ZC> 

PEOUT=PEOUt+1./RR 

VZ2=VZ*VZ 

205 KE=KE+VX*VX+VY*VY*VZ2 

210 CALL PACK(XS (1 ) ,XPACK ( IS) »YPACK (IS) ,ZPACK( IS) ) 

220 PRINT 225»NUM8ER»X*VX*YfVYtZ»VZ»R 
225 FORMAT ( 1H .I0.7E14.7) 

WRITE (jS) XPACK*YPACK,ZPACK 
Q=0 

IF(NS2.EQ.NMK) Q=1 

C DETERMINE WHICH PLOTS IF ANY ARE TO BE MADE. 

GO TO (300 » 240 » 250*260 »270»280* 290) IPLOT 
C MAKE AN X-Y STAR POSITION PLOT. 

240 CALL DDIPLT <0* IN*NBS,XPACK*YPACK*XMIN*XMAX9 YMINtYMAX 
I 3,XPP*1*YP*13» ITAPX) 

GO TO 300 

C MAKE AN X-Z SyAR POSITION PLOT 

250 CALL ODlPLT(Q*IN,NBS*XPACK:,ZPACK,XMIN*XM*Y,ZMIN,ZMAX 
1 3tXPPtlfZP*13*lTAPX) 

GO TO 300 

C- MAKE A Y-Z STAR POSITION PLOT 

260 CALL DDIPLT (Q» IN*NBS*YPACK»ZPACK*YMIN» YMAX»ZMIN»ZMAX 

1 3,YPP»l*ZPil3,lTAPX) 

GO TO 300 

C MAKE A RADIAL VELOCITY VS. RADIUS PLOT 

270 DO 275 IS=1,NBS 

CALL UNPACK { X S ( 1 ) .XPACK ( IS) ,YPACK(IS) *ZPACK(IS) ) 
XC=X-N04 
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YC=Y-N04 

R=SQRt (XC*XC*YC*YC) 

IF (R • EG • 0 • ) GO TO 272 
YPACK < IS ) = { XC*VX+YC*VY > /R 
GO TO 275 

272 YPACK(IS)=:SQRT(VX*VX + VY#VY) 

275 XPACK(lS)=R 

CALL DDIPLT (Q* IN»NBS, XPACK, YPACK, RM IN, RMAX.VRM IN, VRMAX* 

1 3»RPPt 1 *VRP, 13, ITAPX) 

GO TO 300 

C HAKE AN AZIMUTHAL VELOCITY VS. RADIUS PLOT 

280 DO 285 IS=1,NBS 

CALL UNPACK (XS(1) ,XPACK<lS) ,YPACK<(IS) ,ZPACK(IS) ) 

XC=X~N04 

YC=Y-NOA 

R=5SQRT(XC*XC + YC*YC) 

IF ( R • EQ, 0 • ) GO TO 282 
YPACK ( IS) =- (XC*VY-YC*VX) /R 
GO TO 285 
282 YPACK ( IS ) =0 . 

285 XPACK ( I S ) =R 

CALL DDIPLTCQ*IN,NBS, XPACK, YPACK, RMIN,RMAX,VTMIN,VTMAX* 
1 3,RPP,1,VTP»13» ITAPX) 

GO TO 300 

C MAKE AN AXIAL' VELOCITY V5. RADIUS PLOT 

290 CALL SHIFT2(ZPACK(1> ,NBS) 

DO 295 I S = 1 » NBS 
XC=XPACK (IS) -N04 
YC=YPACK(lS)-NOA 
295 XPACK (IS)=SQRT(XC*XC+YC*YC) 

CALL DDIPLT(Q»IN,NBS,XPACK,ZPACK,RMIN,RMAX,VZMIN,VZMAX, 
1 3,RPP»1,VZP*I3»ITAPX) 

■ 300 NS2=NS2+NBS 

NUM9ER=NUMBER+1 
IF(NS2.LT.NBR) GO TO 95 

C REWIND DISK FILES JS AND JT AND EXCHANGE THEIR TAPE NUMBERS 
REWIND JT 
REWIND JS 
jsave=js 


305 


OF CYLINDERS 16) 


C 

C 

c 


TOTAL KINETIC ENERGY* TOTAL 

AND print averaged 


310 


C 

c 

c 

c 

c 


J5=JT 
JT= JS AVE 

FORMAT UH0*40HNUMBER OF particles outside 
IFUPRINT.EQ.I) go to 500 
IF LONG DIAGNOSTIC PRINTING IS TO- B E DONE. PRINT 
POTENTIAL ENERGY AND TOTAL ENERGY. ALSO, COMPUTE 
QUANTITIES as FUNCTIONS OF R* Z» AND R Z. 

DT2=DT*DT 

KE=KE/(8.*DT2) 

PE=(DtE2*PEOUT-PEIN)/DT2 

et=ke+pe 

axial velocity. 

DO 320 I R = 1 . N04 

{R25v6R?IR)=MRTt<VRZAV6H < IIR?-VRAVGR(IR)»VRAVBR<IR>/HASSOR(I«> 

L /(MASSDR(IR)*DT2J) 

St»^1ir7S2X“M 

1 /(MASSDR(IR)*DT2) ) 

/{MASSDR ( IR)*DT2> ) 

VZAVGR ( IR )=VZAVGR(IR)/ (MASSDR ( IR)*DT) 

continue 

MASSDR Cl)=MASSDR<l)*XM/(.25*P I) 

XM02PI = XM/(2*PD 
DO 330 I R“2 » N04 

MASSDR ( IR ) =MASSDR ( IR> *XM02PI/R 
PRINT 340 


320 


330 




*> 


6 £. 




i*v -' : L. * 




'• i-ftt* iiji i Vi.-5.-V t't P- l-t .Ar- « V; 




o o n 


340 FORMAT ( 1H0 *38HRADIUS MASS DENSITY AVG VR* 

1 4BH SIGMA VR AVG VT SIGMA VT* 

2 32H AVG VZ SIGMA VZ/> 

DO 350 IR=1*N04 
R=IR-1 

350 PRINT 360»R»MASSDR<IR)*VRAVGR(IR),VR2AVGR<lR)»VTAVGR(IR)t 
1 VT2AVGR (IR) *VZAVGR ( Ip) ,VZ2AVGR < IR) 

360 FORMAT ( 1H ,F6.1,7E16.8> 

DO 370 KZ=1*NH02 

COMPUTE AND PRINT THE FOLLOWING AVERAGE QUANTITIES AS FUNCTION- 
OF AXIAL POSITION(HEIGHT) - MASS DENSITY* AVERAGE AXIAU VELOCITY* AND 
STANDARD DEVIATION OF AXIAL VELOCITY. 

IF {MASSDZ (KZ) .LT.l.) GO TO 370 

VZ2AVGZ (KZ) =SQRT ( (VZ2AVGZ (KZ) -VZAVGZ(KZ) *VZ A VGZ (KZ) /MASSDZ (KZ) ) 

1 / (MASSDZ (KZ) *DT2) ) 

VZAVGZ(KZ)=VZAVGZ(KZ)/ (MASSDZ <KZ)*DT> 

MASSDZ <KZ)=MASSDZ(KZ)*XM 
370 CONTINUE 
PRINT 380 

380 FORMAT ( 1HO * 42HHE IGHT ( Z ) MASS DENSITY AVG VZ* 

1 16H SIGMA VZ/) 

DO 390 KZ=1 * NH02 

Z=KZ 

390 PRINT 400*Z*MASSDZ(KZ) ,VZAVGZ(KZ) *VZ2AVGZ(KZ) 

400 FORMAT (1H »F6.1*3E16.8) 
c COMPUTE and PRINT AVERAGE RADIAL GRAVITATIONAL! FIELD AS A FUNCTION OF R-Z. 
PRINT 410 

410 F0RMAt(1H0*5QHRADIAL FIELD AS A FUNCTION OF R AND Z (KZ=2*NH02-1 ) ) 

PRINT 420 

420 FORMAT ( 7H RADIUS) 

NH02M1=:NH02-1 
DO 440 IR=2* N04M1 
R~IR-1 

N04P I R=N04+ I R 
DO 430 KZ=2*NH02M1 

430 MASSDZ(KZ)=PHI (N04PIR+I »N04* KZJ-PHI (N04PIR-1 *N04»KZ) 





440 PRINT 445*R* (MASSDZ (KZ) *KZs2»NH02Ml ) 

C 445 COMPU^AND^ AXIAL GRAVITATIONAL 1 FIELD AS A FUNCTION OF R-Z. 

450 FOrUatUHO*45HZ FIELD AS A FUNCTION OF R AND Z (KZ=2.NH02-1 ) ) 

PRINT 460 

460 FORMAT < 7H RADIUS) 

DO 463 KZ=2*NH02M1 

463 MASSDZ (KZ)=PHI (N04*N04,KZ+1 > “PHI (N04*N04*KZ-1 ) 

PRINT 465* (MASSDZ { KZ ) *KZ=2*NH02Ml > 

465 FORMAT (7H 0.0*8E16.8) 

DO 480 IR=2»N04M1 
R=IR-1 

N04PIR=N04+IR 

470 MASSDZ (Kz7=PHl' (N04PIR*N04*KZ + 1 ) -PHI (N04PIR*N04*KZ-1 ) 

480 PRINT 445*R* (MASSDZ (KZ> »KZ=2*NH02Ml ) ,-r H 

C SHIFT THE DENSITY TO THE MOSy SIGNIFICANT HALF OF EACH WORD OF THE PHI MESH. 
500 CALL SHIFT2(PHI(1*I.1) »NCH) 

C WRITE THE DENSITY FROM THE PHI MEPH ONTO DISK FILES 1*2*5* AND t: 

WRITE(l) (({PHI(I*J*K)jI=1*N04)*J=1»N04)»K=1*NH02) 

REWIND 1 

WRITE(H) ( ( (PHI ( I * J*K) * I = N04P 1 » N02 ) *J=1*N04) *K = 1 *NH02) 

WRITE^f ( ( (PHI ( I * J*K) * 1 = 1 *N04) * J=N04P1 *N02) *K=1 »NH02) 

REWIND 5 

WRITE (6) (((PHI(I*U*K) *I=N04P1 *N02) » J=N04P1*N02> *K = 1»NH02> 

REWIND 6 

RETURN 

END 






APPENDIX C 


Computer Plots Produced by the Three-Dimensional Galaxy Simulator 
of Appendix B. 

A long run has yet to be made. These are sample plots showing a 
set of initial conditions and eight subsequent cycles (time steps) . 

The rough appearance of the initial velocity distributions is due to 
the small mesh size (32 x 32 x 8) which necessitated making the galaxy 
only about two cells thick. For an actual run, the dimensions will be 
increased to 64 x 64 x 16. It may prove necessary to generate initial 
conditions in whole or part by analytic means. 


Cycle No. 

Plot Type (Stars in 
Position/Velocity Space) 

Page No 

0 

y position vs x position 

C-2 

0 

z position vs x position 

C-3 

0 

z position vs y position 

C-4 

0 

r velocity vs r position 

C-5 

0 

azim. velocity vs r position 

C-6 

0 

z velocity vs r position 

C-7 

1 

y position vs x position 

C-8 

2 

z position vs x position 

C-9 

3 

z position vs y position 

C-10 

4 

r velocity vs r position 

C-ll 

5 

azim. velocity vs r position 

C-12 

6 

z velocity vs r position 

C-13 

7 

y position vs x position 

C-14 

8 

z position vs x position 

C-1S 
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APPENDIX D 


Listing of the Two-Dimensional Particle-in-Cell Simulator of the Jeans 
Instability in an Infinite Self-Gravitating Compressible Gas. 

Comment cards necessary to make this listing self explanatory will be 
added later. 


Overlay No. 

Program Name 

Page No. 

( 0 , 0 ) 

GASJ 

D-2 

(1,0) 

GETPHI 

D-3 

(2,0) 

INITGAS 

D-6 

(3,0) 

ADVGAS 

D-ll 

(4,0) 

GASPLOT 

D-25 


D 






OVERLAY <IFILE»0»0) 

PROGRAM 6ASj(OUTPUT.TAPEl,TAPE8tTAPE9,TAPE10.TAPEll) 
C0MM0N/ALLC0M/I2 A* ItESt «Ni CY*CYY»RHO (33*33) 

C0MM0N/GA5C0M/NPL0TGtNPRlNTG*'JPRNTMG*N02S0»N02*N02Pl *N02M1 »DMG* 

l GM102* AU02» AU04*SAVPE (200) »SAVIE (200) *SAVKE (200) *SAVTE (200) * 

> CYMlM*N8RG»NBSG»JTG*j5G*NMKG*PI»MAbKl*ING(2 ) *xming*xmaxg* 

} YMING* YMAXG*XPG(2) * YpG* ITAPXG»XYL (3) *PLP (2) *PLT (2) *PLM (2) * 

f RC0S,RSIN,XPR<5) *YPR<5) * XMAX4, XM IN4, YMAX4* YMIN4*SPHIM, SPHlT t 

> S PH IP, EMIN, EMAX,ScALM,DT*DT2,DT3,DT4,SUMMAS, SUMP*, SUMPY* 

b Sl)MKE*SUMIE,SPHl*PL (2) 

INTEGER CY * CYY 

IFILE-5LIFILE 

GPFILE=6LGPFILE 


AGFILE=6LAGFILE 

CALL OVERLAY ( IFILE*2»0*6NRECALL) 
ITE5Tsl 

10 PRINT 5 * CYY 
5 FORMAT ( 1H0»6HCYCLE = * 18) 

CALL SECOND (T1 ) 

CALL OVERLAY (GPFIlE, 1 *0*6HRECALL) 
CALL SECOND ( T2 ) 

T3=T2-Tl 
PRINT 15* T3 

15 FORMAT ( 12H FIELD TIME=,E16.B> 

• CALL OVERLAY (AGFILE*3,0*6HRECALL) 
CALL 5EC0ND(T4) 

T5=T 4-T2 
PRINT 35*T5 

35 FORMAT (14H ADVANCE TlME=, E16.S) 
IF(CYY.GE.CY) GO TO 40 
CYY=CYY+1 
GO TO 10 

40 CALL OVERLAY (IFILE*4»0 »6hRECALL) 
END 


\ 

N 


OVERLAY (GPFILE* 1 »0) 
program GETPHI 

COMMON/ ALLCOM/I2A,IlESTtN*CY»CYY*RH0{33t 33) 
COMMON Z (1025) *Y (1025) 

DIMENSION G ( 33 » 33) 

N02-N/2 

M21=N02+1 

IF(ITEST.EQ.O) GO TO 40 
I 2B= I 2A- 1 
RNI=1./(N*N> 

N02P1=N21 
DO 1 J~ 1 *N 
DO 1 1 = 1 *N 

IF(I.EQ.l.AND.J.EO.l) GO TO 1 
G(I,J)=RNI/SQRT< (I-1.)*CI-1.)+(J-1*>*(J-1.> ) 

1 CONTINUE 
G<l*l)-G(l*2) 

CALL GETSET (2» I2B) 

DO 3 J=1,N02P1 

DO 4 1=1 »N02P1 
4 Z ( I ) =G ( I » J ) 

CALL FTRANS(2,I2B) 

DO 10 1=1 *N02P1 
10 G ( I * J ) =Y ( I ) 

3 CONTINUE 

DO 2 1=1 »N02P1 
DO 19 J=l t N02Pl 
19 Z(J)=G(I.J) 

CALL FTRANS (2* I2B) 

DO 22 J=1 *N02P1 
22 6(ltJ)sY(JJ 

2 CONTINUE 
WRITE(l) G 
REWIND 1 

40 READ ( 1 ) G 
REWIND 1 

CALL GETSET (3» I2A) 

DO 7 J=ltN 
DO 8 1 = 1 *N 
8 Z<I)=RH0<ItJ) 


CALL FTRANS ( 3 * I2A) 

DO 9 1=1, N 

9 RHO(I,J)=Y(I) 

7 CONTINUE 
DO 25 1=1, N 
DO 26 J=1,N 

26 Z< J)sRHO(ItJ) 

CALL FTRANS (3, I2A) 

DO 27 J=1*N 

27 RHO < I , J) =Y < J) 

25 CONTINUE 

DO 11 J=2 , N02 
DO 11 I=2»N02 
I1=I + N02 
J1=J*N02 

RH0(I,JJ=RH0(I*J)*6«I,J) 
RHO(Il,J)=RHOUl, J>*G<I,J) 

RHO ( I,Jl)=RHO(I, J1 )*G< I,J) 

RHO ( 1 1 ,U1 ) =RHQ ( 1 1 * J1 ) «G ( I , J > 

11 CONTINUE 

DO 37 1=2, N02 

RHO ( I ♦N02 ,1)=RH0(I + N02,1)*G(I,1) 

RHO ( I +N02 ,N21 ) =RHO ( I +N02, N21 ) *G ( I »N21 ) 
RHO(I,l)=RHO(I,l)*G<I,l) 

37 RHO(I,N21)=RHO<I,N21>*G(I,N21) 

DO 38 J=2,N02 

RHO (1 , J ) =RHO ( 1 » J ) *G ( 1 , J ) 

RHO ( 1 , J+N02 ) =RHO ( 1 » J+N02 ) *G ( 1 » J) 

RHO { N21 * J) =RHO (N21 , J) *G(N21 , J) 

38 RH0(N21, J+N02)=RH0(N21, J+N02) »G<N2l , J) 
RHO ( 1 , 1 ) =RHO (l,l) tt G(l,l) 

RHO (N21 , 1 ) =RHO (N21 , 1 ) *G <N2 1 » 1 ) 
RH0(1,N21)=RH0(1,N21>*G<1»N21) 
RH0(N21,N21 J=RH0(N21,N21 )*G(N21,N21) 
CALL GETSET (A, I2A) 

DO 14 J=1,N 
DO 12 1=1, N 

12 Z(I)=RHO(I,J) 

CALL FTRANS (A, I2A) 

DO 13 1=1, N 


13 RHOU*J)=Ym 

14 CONTINUE 

DO 2B 1=1 *N 
DO 29 J=1*N 

29 Z ( J) =RHO ( I ♦ J) 

CALL FTRANS(4,I2A) 
DO 30 J=1*N 

30 RHQ ( I » J) -Y ( J) 

28 CONTINUE 
DO 81 1 = 1 * N 
RHO(N*l*I>=RHOU»I> 
81 RHO<I*N*l)=RHO(I*l> 
ITEST=0 
RETURN 
END 


ru n 


SET 


OVERLAY (IF1LE»2*0) 

PROGRAM INITGAS 

g83HO^feSS^5KfNJSiars:5?«^!JS:3S5iSiMoe.NO*P». S oa ; . T o ; .6. 

1 GM 102 ♦ AU02 * AUOAt SftyPE ( 200 ) *SAVlE(200) *SAVKE(200) *5AVTE(200) * 

CYMlNtNB'lGtNBSG* JTGt JS»G«NMKG*PI*MASK1* IN6(2) *XHING*XMAXGt 

YMING* YMAXGtXPG (2) *YPG, ITAPXG*XYL (3) »PLP<2) »PLT(2) » 

4 RC0S,RSIN,XPR(5).YPR{5) ,XMAX4* XM IN4* YMAXA* YM IN4, SPHlM* SPHlT . 

5 SpHlP«EMINtEMAX*SCALM*DT*DT2*DT3«DT4»SUMMAS t SUMPX*SUMPY* 

6 SUMKE*SUMIE*SPHI»PL(2) Urtlf w C oa % 

DIMENSION A 1 (32*32) »A2 (32*32) * A3 (32*32) *A4 (32*32) *XH0L0(528> ♦ 

^EQUIVALENCE^^RHO ( 1 * 1 ) * XHOLD ( 1 ) ) t (RHO (1*17) *YHOLO( 1 ) ) 

INTEGER CYtCYY 

REAL IN I T IE „ 

INITIAL VALUES 

SET SIZE WHERE N=N02-2*#I2A IS THE SIZE OF THE ACTIVE MESH 
I2A=5 

SET TOTAL NUMBER OF GAS PARTICLES 

NBRG=1 2672 „ _ 

Set number of gas particles per read of gas PARTICLE FILE 

NBSG=528 

SET ARTIFICIAL VISCOSITY CONSTANT OF GAS 


C 


C 

C 


C 


C 

C 





C 

C 


AU=Q. 

SET MASS PER GAS PARTICLE 
XMG=1 . 

SET THE INITIAL SPECIFIC INTERNAL! ENERGY 
INITIE=5. 

SET TOTAL NUMBER OF TIME 5tEP5 
CY=150 

SET NUMBER OF POINTS PER GAS PLOT 


npg=nbrg 

SET PERIOD OF GAS LONG PRINTING 
NPRINtG=50 

SET PERIOD OF GAS DENSITY SHORT PRINTING 


NPRNTMG-25 

SET PERIOD OF GAS PLOTTING 


5ET°RATI0 OF SPECIFIC HEAT AT CONSTANT VOLUME TO 
PRESSURE. (GAMMA MUST BE GREATER THAN OR EQUAL! TO 


SPECIFIC HEAT AT CONSTANT 
1. GAMMA IS EQUAL- TO 2.0 


o o 


FOR NORMAL SIMULATION OF MONATOMIC GAS IN TWO DIMENSIONS. GAMMA MAY BE SET 
EQUAL TO 1.0 FOR A SIMULATION WITHOUT PRESSURE TERMS.) 

c set M itape=i to start run. set itape =2 to continue run with pick up tape 

ITAPE=1 

C SET CONSTANTS 

PI=3. 1415926536 

MASK 1=0777777777700000 000 0 0 

JTG=9 

JSG-10 

GM 102= ( GAMMA-1 •) /2. 

AU02=AU/2. 

AU04=AU/4. 

NMKG=NPG“NBSQ 

N=2**I2A 

c note following definition 

N02=N 

N02P 1 =N02* 1 
NQ2M1 =N02-1 
N02SQ=N02*N02 
N04=N02/2 

GMU=XMG*NBRG/N02SQ 

tau=1/Sqrt (gmu> 

DT=T AU/50 • 

DT2=DT*DT 

DT3=DT2*DT 

DT4-=DT2«DT2 

C SET PLOTTING SPECIFICATIONS 
ING<1)=10H2D-GAS 
ING ( 2 ) =1 OH 
XMING=-10 
XMAXG=40 
YMING=- 1 0 
YMAXG=40 

XPG ( 1 ) =1 OHX * CYCLE- 

YPG=1 OH Y 

ITAPXG=6LTAPE23 

XYL ( 1 ) =1 OHX-Y PLANE, 

XYL(2)=10H CYCLE = 

PL(1)=10HPOTENTIAL* 

PLP ( 1 ) =1 OH PRESSURE* 


PLT < 1 > =7. OH 1EMP* 

PLM ( 1 ) s l OH DENSITY* 
ANG=Pl/3. 

RC0S=C0S (ANG) 

RSIN=SIN (ANG) 

IC=N02+1 

XPRU)=1.+N02*RC0S 
XPR ( 2) =1 • 

XPR (3 ) si * *N02 
XPRU) = IC + N02*RC0S 
XPR ( 5 1 =XPR { 1 ) 

YPR ( l ) =1 . ♦NOSERS IN 
YPR ( 2 ) = 1 . 

YPR ( 3 ) - 1 • 

YPR (4) =1 •♦N02*RSIN 
YPR (5) =YPR ( 1 ) 
XMAXA=IC+N02*RC0S 
XM INA-=- 1 . 

YM IN4=- I . 

YMAX4=XMAX4 

HT=IC+N02»RC0S-N02*RSIN 
■ PMAX=2.*NBRG*XMG*DT2/N04 
PHAXM=10.*GMU*DT2 
PMAXT=50.*INITIE*DT2 
PMAXP=.2*PHAXM*PMAXT 
SPHI=HT/PMAX 

ENCODE { 1 0 » 25»PL ( 2) > SpHI 
SPHI M=HT/PMAXM 
ENCODE ( 1 0 *25* PLM (2) ) SPHIM 
25 FORMAT (F10. 3) 

sphit=ht/pmaxt 

ENCODE (10*25»PLT(2)> SpHIT 
SPHIP=HT/PMAXP 
ENCODE <10*30* PLP(2) ) SPHIP 
30 FORMAT (FIO.O) 

EMIN-- ( <XMG*NBRG) **2) /NOA 

EMAX=-EMIN 

SCALM=GMU*DT2/7. 

INITIE=INITIE*DT2 

DMG=XMG*Dt2 

DMG=DMG* AND .MASK 1 


IF ( ITAPE.EQ.2) GO TO 350 
CYY-1 

DO 105 1=1. N02 
DO 105 J=1 . N02 

AM I. J)=0. 

105 CONTINUE 

X=URAN (7654321 « ) 

NS2=0 

110 DO 130 IS=1,NBSG 
X=N02*URAN ( 0 • ) ♦ .5 
Y=N02*URAN (0 • ) ♦ *5 
I X=X + . 5 
JY=Y ♦ .5 

IFUX.GE.l) GO TO 115 
X=X+N02 
IX = X + ,5 
GO TO 120 

115 IFIIX.LE.N02) GO TO 120 
X=X-N02 
IX = X + • 5 

120 IF ( JY.GE. 1 ) GO TO 125 
Y = Y *N02 
JY=Y+ *5 
GO TO 128 

125 IFUY.LE.N02) GO TO 128 
Y«Y-N02 
JY=Y*.5 

? 28 AMIX.JY)=AMIX.JY)*DMG 
XHOLD ( IS) =X 
YHOLDUS)=Y 
130 CONTINUE 

WRITE (9) XHOLD. YHOLD 
NS2=N52*NBSg 

IF(N52.LT.NBRG) GO TO 110 

REWIND 9 

SUMMAS=0 

5UHIE=:DMG*NBRG*INITIE 

5UMKE=0. 

surpx=o. 

SUMPY=0. 

\ 



DO 140 1=1 *N02 
DO 140 J=ltN02 
IF(A4<I»J) »EQ» 0 • ) GO TO 135 
5UMMAS=SUMMAS+A4( I»J) 

A3 < I • J)=INlTlE 
GO TO 130 
135 A3 { I * J ) =0 
138 A1 ( I ♦ J) =0 • 

A2 ( I » J ) =0 • 

RH0(I*J)».5*A4(I*J> 

140 CONTINUE 

DO 150 1=1 *N02P1 
RHO ( I tN02Pl ) =0 * 

RHO (N02P1 * I ) =0 • 

150 CONTINUE 
GO TO 400 

C STATEMENTS 350 TO 400 ENABLE RUN TO BE CONTINUED WITH PICK- UP TAPEo 
350 NS2=0 

355 READ (11) XHOLD*YHOLD 
WRITE (9) XHOLDt YHOLD 
NS2=NS2+NBSG 

IF(NS2.LT.NBRG) GO TO 355 
REWIND 9 

READ ( 11 ) A1*A2*A3*A4*CYY* SUMMAS * SUMP X* SUMP Y »SUMKE*SUMIE 

CY=CY+CYY 

CYY=CYY ♦ 1 

DO 360 1=1 *N02 

DO 360 J=1*N02 

RH0(I»J)=.5*A4(I*J> 

360 CONTINUE 

DO 370 1=1 *N02P 1 
RHO(I*N02P1)=0. 

RHO ( N02P 1 * I ) =0 • 

370 CONTINUE 
400 CYMIN=CYY-1 

CALL EVICT (6LTAPE11) 

C WRITE U»V*I* AND M ONTO TAPES 
WRITE18) AlfA2.A3»A4 
REWIND 8 
RETURN 
END 


OVERLAY (AGFILE*3*0> 

PROGRAM ADVGAS 

COMMON/ALL COM/ 1 2A* UESy *N*CY*CYY*RHO (33*33) 

COMMON/GASCOM/NPLOTG* NPR INTG* sjPRNTMG*N02SQ ,N02*N02P 1 « N02M1 *DMG* 

1 GM102 * AU02 » AU04 *SAVPE ( 200 ) ,SAVIE(200> * SAVKE ( 200 > ♦ SAVTE ( 200 ) , 

2 CYMIN.NBRG*NBSG»JTG*JSG»NMKG»PI*MASK1 * ING ( 2 ) *XMlNG*XMAXG* 

3 YM ING* YMAXG*XP6 (2 ) *YPG*lTAPXG*XYL(3) *PLP(2> *PLT (2) *PLM (2) * 

4 RC0S,RSIN*XPR<5) * YPR ( 5 ) *XMAX4, XM IN4* YMAX4* YM IN4* SPHIM * SPHIT * 

5 SPHIP*EMIN*EMAX»SCALM*DT»DT2»DT3»DT4*SUMMA5,SUMPX*SUMPY, 

6 SUMKE*SUMIE*SPHI»PU2) 

DIMENSION A1 (32*32> *A2 (32*32) * A3 (32 *32) *A4 (32.32) ,AlH0R2 (32) * 


1 A1VERK32) *A2H0R2(32) *A2VERT<32) *A3HORZ(32) *A3vERT(32) . 

2 XHOLD (528) .YHOLD (528) *PX ( 32* 32) *PY (32*32) *ET0TAL (32*32) 


EQUIVALENCE ( RHO ( I * l ) ♦ XHOLD ( 1 ) ) ♦ (RH0(1*17) ,YHOLD(l) ) 
INTEGER CYY *CY*Q 

REAL MC*MIP*MIM*MJP,MJM. IE*KE, IC 


REAL IEIP*IEIPJP* IEIPJM* IEJP* IEJM*IEIM»IEIMJP*IEIMJM 
C If GAS PLOTTING IS JO BE DONE THIS CYCLE SET IPL0TG=1» OTHERWISE SET IPL0TG=Q 
IPLOTG-O 


IF (CYY-CYY/NPLOTG»NPLOTG.EQ.O.OR.CYY.EQ*1 ) IPLOTG=l 


CALL PSEUDO 

ENCODE (10*21* XPG(2) ) CYY 
ENCODE (10*21*XYL(3) ) CYY 
21 FORMAT(IIO) 

IF< IPLOTG.NE.l) GO TO 120 

C MAKE A CONTOUR PLOT OF THE GRAVITATIONAL POTENTIAL' 

CALL DOIPLT<0*ING*5*XPR*YPR*XMIN4*XMAX4*YMIN4*YMAX4*3*XYLU ) * 
1 2*PL(1)*14.ITAPXG) 

DX=0. 

DY = 1. 

DO 50 J=2*N02 


K=0 

DX=DX*RC05 

DY=DY*RSIN 


DU=0. 


DV-1 . 

DO 40 I=2*N02 
K=K + 1 

DU=DU+RCOS 

dv=dv+rsin 


A1H0RZ(K)=J+DU 
AIVERT(K) =SPHI*RHO( J» I) *DV 
A2H0RZ (K > = 1 + DX 

40 A2vERT(K>=SPHI»RHO(I*J) *DY 

CALL DDIPLT <0* ING* K » A2H0RZ* A2 VERT « XMI N4* XMAX4. YM I N4* YMAX4 t 

1 3*XYL(l)»2fPL(l)»14»ITAPXG) 

CALL DDIPLT (0» ING»K»A1H0RZ*A1VERT*XMIN4*XMAX4*YMIN4*YMAX4* 

1 3*XYL(1)»2*PL(1) »14-.ITAPXG) 

50 CONTINUE 
K=0 

DO 60 I=2*N02 
XI = I 
K=K ♦ 1 

A2H0RZ (K)=XI+RCOS 
A2VERT<K»=1.+RSIN 
60 CONTINUE 

CALL DDIPLT ( 1 » ING»K * A2H0RZ * A2VERT * XMIN4*XMAXA* YM IN4» YMAX4* 

1 3,XYL(1)*2*PL(1)*13*ITAPXG) 

C IF GAS PRINTING IS TO BE DONE THIS CYCLE SET IPRlNTG=l tOTHER^lSE SET IPRINTG-0 
120 IPRINTG=0 

IF(CYY-CYY/NPRINTG*NPRINtG.EQ.O.OR.CYY.EQ. 1 ) IPR INTG=1 
C IF GAS DENSITY SHORT PRINTING IS TO BE DONE THIS CYCLE SET IPRNTMG=1* 

C OTHERWISE SET IPRNTMG=0. 

IPRNTMG-0 

IF {CYY-CYY/NPRNTMG*NPRNTMG.EQ.0.0R.CYY.EQ.1 ) IPRNTMG-l 
C READ U»V* I »AND M FROM TAPE8 INTO Al»A2.A3. AND A4 RESPECTIVELY 

READ (8) Al *A2* A3» A4 
REWIND 8 
nhivel=o 
sumpe=o. 

DO 200 J=1*N02 
DO 200 I “1 *N02 
MC~A4 { I * J> 

IFIMC.LT. DMG) GO TO 183 
5UMPE=SUMPE+MC*RH0 ( I ♦ J) 

UC = A 1 { I * J > 

VC=A2 { I* J) 

PC=A3( I, J)*MC 

IF ( I • EQ. 1 ) GO TO 130 
IF ( I.EQ.N02) GO TO 140 



i'/-y 


UIP=A1 ( 1 + 1 * J) 

MIP=A4(I+l*J) 
PIP=A3(I+l*J)*A4(I*l*J) 
UIM=A1 (I-l* J> 

M IM-A4 ( I - 1 * J > 
PIM=A3(I-1*J)*A4<I-1*J) 
GX=RHO (I*1*J)-RH0(I-1*J) 

GO TO 150 

130 UIP=A1(2,J) 

MIP=A4(2»J) 

?1P=A3(2*J)*A4(2*J) 

UIM=A1 (N02» J) 

HIM=A4(N02» J) 
PlM“A3(N02*J) tt A4(N02*J) 
GX-RHO ( 2 * J) -RHO (N02* J) 

GO TO 150 
140 UIP=Altl,J) 

MIP=A4(1*J> 

PIP=A3(1*J)*A4(1.J> 

UIM=A1 (N02M1 * J) 

MIM=A4 (NQ2M1 ♦ J) 

PIM=A3(N02M1*J)*A4(N02H1*J) 

GX=RH0(l*J)-RKO?N02MUJ) 

150 IF(J.EQ.l) GO TO 160 
IF(J.EQ.N02) GO TO 170 

VJP-A2 { I * J*1 > 
mjp=A4(I»j+i:< 
PJP=A3(I*J*1)#A4(I*J*1) 
VJM=A2(I,J-1 ) 

MJM=A4 ( I * J-l ) 
PJM=A3U*J-l)*A4(If J-l) 

GY=RHO (I ♦ J* 1 ) -RHO ( I * J-l ) 

GO TO 175 
160 VJP=A2 (1*2) 

MJP=A4 (1*2) 

PJP-A3(I*2)*A4(I*2) 

VJM=A2(I*N02) 

MJM=A4(I.N02) 

PJM-A3 ( I »N02 > *A4 ( I *N02) 
GY=RH0(I*2)-RH0(I*N02> 

GO TO 175 



170 VJP=A2(Ifl) 

MJP=A4(I»1> 

PJP=A3U»l)*A4(Itl> 

VJH=A2 ( I *N02M1 ) 

MJM=A4 ( 1 *NQ2M1 > 

PJM=A3 ( I »N02M1 ) *A4 ( I »N02Ml ) 

GY=RHO ( I * 1) -RHO < I ,N02Ml ) 

175 QIP-MIP* (UC-UIP) / (MOM IP) 

IFCQIP.LE.O. ) QIP=0. 

3IM=MIM*(UIH-UC)/(MC+MIM) 

IFtQlM.LE.O.) QIM=0« 

GM=GM1 02/MC 
IF(PC,LT.0.) PC=0. 

PlPB=POPlP 

PlMB=POPIM 

PJPB=PC+PJP 

PJM6=POPJM 

IFtPlP.LE.O.) P IPB=0 • 

IFCPIM.LE.O.) PIMB=0. 

IF (PJP«LE«0* ) PJPB=0. 

IF ( P JM .LE . 0 . > PJMB=0. 

UT=UC-»GH* ( P I MB-P IPB) +AU02*(Q1M^QIP) +GX 
QJP-MJP* (VC-VJP) / (MC+HJP) 

IF {QJP.LE.O • > QJP=0 • 

QJM-MJM* (VJM-VC) / (MC+MJM) 

IF (Q JM «LE • 0 « ) QJM = 0 • 

VT=VC+GM* (PJMB-PJPB) +AU02* (QJM-QJP) *GY 

UIPB-UC+UIP 

UIMB=IK>UIM 

VJPB-VC+V JP 

VJMB=VC*VJM 

IEsA3( ItJ) ♦ .5*(UC*UOVC*VC* (UC+UT ) *GX» ( VC+VT > *GY ) 

1 ♦AUQ4*<QIM*UIMB-QIP*UIPB*QJM»VJM8-QJP»VJPB) 

2 +GM*(PIMP*UIMB-PlP8*UIPB+PjMB*VJMB-PJPB*VJPB)/2. 

3 -.5*(liT*UT + VT ft VT) 

IE=tfC*IE 

IF((^C,LT.5.0*DMG) 60 TO 1B0 
IF(UT.GE..375.0R.VT.6t..375) NHI VEL=NHIVEL+1 
180 IF (ABS (UT> »LT* 1 .0 .AND. ABS { VT) *LT»1 .0) GO TO 185 
. CY=CYY 



IPL0TG=1 

IPRINTG=1 

JPRMTMG=1 

PRINT 182 »Ut *VT » I y J 

182 FORMAT ( 1H0*UT=*E16«8* AND VT=*E16.8* FOR I=*I2* AND J=*l2) 
GO TO 185 

183 UT=99. 

VT=99, 

IE=0 • 

185 PX(I»J)=UT 


200 


205 


PY ( I * J) -VT 

ETOTAL(I.J)=IE 

CONTINUE 

IF(NHIVEL..GT.0) PRINT 205*NHIVEL 

Format (1H0*NUMBER of cells containing greater than 5 particles and 
1 HAVING SCALED VELOCITY COMPONENT (S) GREATER THAN ,375= *15) 

DO 210 I=ltN02 


DO 210 J-l *N02 
A1(I#J>=PXU*J) 


PX(I» J)=0. 

A2 ( I * J ) =P Y ( I f J ) 
PY<I*J)=0. 

A3 ( 1 1 J) ”ETOT AL d » J) 

ETOTALdf J)=0. 

210 CONTINUE 


DO 222 I“1 yN02 


DO 222 J=l»N02 
222 RHO ( I i J) -0 • 

DO 224 I=2*N02Ml 
DO 224 J=2.N02M1 
IF ( A3 ( I « J) .GE. 0 . ) GO TO 224 

IE-A3 ( I y J) 

IEIP=A3<I+ly J) 

IF { IEIP.LT *0 • ) IEIP=0. 
IEIPJP=A3<I+Iyj+1) 
IFUEIPJP.LT.O.) IEIPJP=0. 
tEIPJM=A3ddyJ-l> 
IF(IEIPJM.LT*0.) IEIPJM-0 • 
IEJP-A3 < 1 1 J*1 ) 
IFdEjP.LT.O.) IEJP=0* 


IEJM-A3 M* J-l > 

IF ( IEJM*LT«Q * ) IEJM=0. 

IEIM=A3(I-1»J) 

IF { IEIM *LT • 0 • ) IE IM=0 , 

IEIM JP- A3 < I - 1 ♦ J ♦ 1 ) 

IF (IEIMJP.LT *0 • ) IEIMJP-0. 

IEIMJM-A3M-1 ♦ J-l ) 

IFdETMJM.LT.O.) IEIMJM=0. 

sum-Ieip*ieipjp+ieipjm+iejp+iejm+ieim+ieimjp*ieimjm 

IF ( SUM • LT • - 1 E ) GO TO 223 
QUOTNT=:IE/SUM 

RHO(I+l*J)=RHO(I*l » J) ♦QUoTNT*IEIP 
RHOM + lfJ + DsRHOd + lfJ*!) +QUOTNT* IEIPJP 

RHO ( I *1 * J-l ) s RHO < !♦ 1 * J-l ) ♦QOOTNT* IEIPJM 
RHO ( 1 1 J + l ) -RHO ( I« J+i ? ♦QUOTNT* IEJP: 

RHO < I • J-l > =RHO ( I * J-l i ♦QUOTNT* IE JM 
RHO < 1-1 ♦ J) =RHO ( 1-1 * J) +QUOTNT* IEIM 1 

RHO ( 1-1 f J*1 ) =RHO ( 1-1 * J*1 ) +QUOTNT*IEIMJP 
RHO <1-1* J-l) =RHO ( 1-1 * J-l J ♦QUQTNT* IEIMJM 

223 RHO { I * J) =RHO ( I * J) -IE 

224 CONTINUE 

DO 228 X=1*N02 
DO 228 J=1 * N02 
IF (A4 ( I ♦ J) •LT-DMG) GO TO 226 
A3<I*J)=<A3(i*J)*RH0<I* J) )/A4(I*J) 

IF{A3(I*J> .LT.O.) A3(I.J)=0. 

GO TO 227 

226 A3 { I » J ) =G * 

227 A4CI»J)=0. 

228 CONTINUE 
PXX=SUMPX/DT3 
PYY=SUMPY/DT3 

' TM=SUMMAS/DT2 
PE=-SUMPE/DT4 
IE=SU H ie/DT4 

K£=SUMKE/DT4 
TE=PE*IE*KE 
PRINT 232 

232 FORMAT (1H0 *THE VALUES ON THE NEXT TWO LINES ARE FOR THE END OF TH 
IE LASt CYCLE*) 


TM=*E14.7) 


234 

236 


250 


PRINT 234fPXX»PYY,TM 

FORMAT ( 1H *PX=«E14.7* PY=*E14.7* 


PRINT 236tPEtIE*KE*TE 
FORMAT ( 1H ,*PE=*E14.7* 
ICYPUT=CYY-CYMIN 


IE=*E14.7* 


KE-*E14»7° 


SAVPE (ICYPLT)=PE 
5AvIE(ICYPLT>=IE 
SAVKE(ICYPLT)-KE 
5AvTEUCYPLT)=TE 
PRINT 250 

FORMAT (lH0t47HLAST PARTICLE OF EACH SET OF NBSG GAS 


PRINT 260 

260 FORMAT OH ^NUMBER X 

1 VEFF*> 

NS2=0 

NUMBER=1 

270 READ(JTG) XHOLDtYHOLD 
DO 350 IS=1,NBSG 
X=XHOLD( IS) 

Y=YHOLO(IS) 

IX0L0=X+.5 
JYOLD-Y + »5 


IX=X 

JY=Y 

DX-X-IX 

DY=Y-JY 

D11=1-DY-DX+DX*DY 

D12=DYM1-DX> 

D2 1-DX* ( 1-DY ) 


D22=DX*DY 

IFUX.LT.l) IX=IX+N02 
IF(JY«LT.1> JY=JY +N02 

IXPI=IX*1 

JYPI=JY*1 

IFIIXP1.GT.N02) IXP1=IXP1-N02 
IF { JYP1 «GT *N02> JYP1=JYP1-N02 
IF(A1 (IX*JY) .GT.98.) D11=0. 

IF ( A1 (IX*JYPX) .GT.98.) 012=0. 
IF(A1 (IXPltJY) .GT.9B.) 021=0. 


TE=*E14.7) 


PARTICLES) 

UEFF 


DSUM=D1 1 + D 12 +D2 1+022 

UEFF=(Dll tt Al ( IX» JY) +D12*Al (IX* JYP1) +D21*A1 (IXP1*JY> 
1 ♦D?2*A1 ( IXP1 * JYP1 ) ) /D^UM 

VEFF= ( Dl 1*A2 ( I X * JY) +D12*A2 ( IX» JYP'l ) +D21*A2 « IXP1 * JY) 
1 +D22*A2(IXP1*JYP1) J/DSUM 


X-X + J EFF 
Y=Y+VEFF 


lXNEW=X+.5 

JYNEW=Y+.5 

IF ( IXNEW„GE.l) GO TO 310 


X=X + N02* 

IXNEW=X*.5 
GO TO 320 

310 IFUXNEW.LEoN02) GO TO 320 

X=X-M02 

IXNEW=X*.5 

320 IF(JYNEW.GE.l) GO TO 330 

Y=Y*M02 
JYNEW=Y+ .5 
GO TO 340 

330 IF(JYNEW.LE.N02> GO TO 340 
Y=Y-N02 


340 


350 


360 

365 


JYNEW = Y + *5 

DPX=DMG*A1 ( IXOLD* JYOLD) 

DPY=DMG#A2 (IXOLD* JYOLD) 

DE=DMG*< A3 ( I XOLD* JYOLD) *.5* ( Al < IXOLD* JYOLO) *A1 ( IXOLD* JYOLD) 

1 +A2 (IXOLD* JYOLD) *A2< IXOLD, JYOLD) )) 

PX (IXNEW* JYNEW) -PX( IXNEW, JYNEW >*DPX 
PY( IXNEW, JYNEW )=PY (IXNEW, JYNEW) +DPY 
ETOTAL(lXNEW,JYNEW)*ETOTAL (IXNEW, JYNEW) +DE 

A4 (IXNEW* JYNEW )=A4( IXNEW, JYNEW) + DMG 
XHOLD { 1 5 ) =X 
YHOLD ( IS ) sY 

WRITE(JSG) XHOLD* YHOLD 
IF(NUMBER.GT.IO) GO TO 365 
PRINT 360 *NUMSER» X* Y ♦ UEFF * VEFF 
FORMAT ( 1H *I6,4E16.8> 

NUMBER=NUMBER+1 
IF(IPLOTG.EQ.O) GO TO 370 
Q=Q 


IF(NS2.EQ.NMKG) Q=1 

CALL DDlPLT{Q,IMG,NBSG,XHOLD»YHOLO»XMINGtXMAXG»YMING*YMAXG f 
1 2*XPG*l*YPGfl3»lTAPXG) 

370 NS2=NS2+NBSG 

IF{N52.LT.NBRG) GO TO 270 
REWIND JTG 
REWIND JSG 

C EXCHANGE TAPE NUMBERS OF JTG AND JSG 
JSAVE=JSG 

jsg=jtg 

jtg=jsave 

C COMPUTE .5*DENSITY AND STORE IN RHO 

C COMPUTE UP AND STORE IN Al. COMPUTE VP AND STORE IN A2. COMPUTE IP AND STORE 
C IN A3. 

SUMMA5=0. 

SUMPX=0. 

SUMPY=0. 

SUMKE=0. 

sumie=o. 

DO 380 1=1 »NO? 

DO 380 J=1»N02 
MC=AA( I, J) 

IF(MC.EQ.O.) GO TO 376 
RHO ( I * J) = .5*MC 
SUMMASsSUMMAS+MC 

SUMPX-SUMPX+PX( I» J) 

SUMPY=:SUMPY + PY{I» J) 

Al (I»J)=PX(ItJ)/MC 
A2(I»J)=PY(ItJ)/MC 

KE=.5*MC*(A1 (I* +A2 ( I » J) tt A2 ( I » J) ) 

sumke=sumke+ke 

IE=ETOTAL ( It J) -KE 
SUMIE=SUMIE*IE 
A3(I»J)=IE/MC 
GO TO 380 
376 A1(I,J)=0. 

A2 ( I » J ) =0 • 

A3 ( I t J) =0 • 

RHOUt J)=0. 

380 CONTINUE 

DO 390 I=lfN02Pl 


RHO ( I ,N02P1 ) =0 • 

RHO ( N02P 1 * I ) =0 . 

390 CONTINUE 

WRITE18) A1 , A2 » A3, A4 
REWIND 8 

C IF GAS PRINTING IS NOT TO BE DONE THIS CYCLE* GO TO 485 
IF(IPRINTG.NE.l) GO TO 485 
400 FORMAT ( 1 H ,I3*8E15.7) 

N02M7=N02-7 • 

DO 420 IMIN-1 *N02M7 *8 
IMAX=IMIN*7 

PRINT 4IO*IMIN,IMAX*N02 

410 FORMATI1HO*J ( (U< I* J) * I=*I3*,*I3*> * J=1,*I3*>*) 

420 PRINT 400* (J* (A1 (I»J) *I=IMIN.IMAX)*J=1*N02) 

DO 440 IMINsl ,N02M7,8 
IMAX=IMlN+7 

PRINT 430* IMIN* IMAX»N02 

430 FORMAt( 1HO»J < < v < I * J) * I=*I3»»*I3« ) » J»1 *«I3*) «) 

440 PRINT 400* (J* <A2(I*J) *I = IMIN*IMAX) ,J=1»N02) 

DO 460 I M I N=1 * N02M7 * 8 
IMAX= JM IN+7 

PRINT 450, IMIN. IMAX.N02 

450 FORMAt(1H0*J ( ( I C I * J) * I»*I3***I3* > * J=X»*I3*) *> 

460 PRINT 400, ( J, (A3 ( I , J) ♦ I=IMIN, IMAX) »J=1,N02) 

DO 480 IM1N=1,N02M7,8 
IMAX=IM IN+7 

PRINT 470, IMIN, IMAX, N02 

470 F0RMAT<1H0*J ( (M ( I , J ) , I=* 13* ,* 13* ) ♦ J=1*I3* ) * ) 

480 PRINT 400, (J, (A4(I,J) ,I=IMIN,IMAX)*J=1.N02) 

C IF GAS PLfTTlNG IS NOT TO BE DONE THIS CYCLE, RETURN. 

485 IF(IPLOTG.NE.l) GO TO 600 
C PLOT VELOCITY VECTORS 
XN02=N02 
A = 0 • 

CALL DDIPLT(0, ING, 1 , A, A* 0 , XN02 , 0 ,XN02,2,XPG, 1 , YPG, 13, ITAPXG) 
XYSCALE=10./NO2 
VSCALE=. 1/DT 
' DO 500 1=1, N02 
DO 500 J=1 *N02 

IF(A4(I,J) .LT.DMG) 0 ^ TO 500 


i 




n> 

I ! XA=XYSCALE*I 
YA=XYSCALE*J 

XB=XYSCALE»( I+VSCALE»A1 (I» J) > 

YB=XYSCALE*< J+VSCALE»A2(I* J) ) 

CALL PARR0W(XA,YA»XB*YB»1> 

500 CONTINUE 

CALL 0DIPLT<1*ING,1 tA,A*0»XN02tO»XN02»2»XPG.l*YPG«13, ITAPXG) 

C MAKE A CONTOUR PLOT OF THE DENSITY 

CALL DDIPLT<0,ING,5*XPR,YPR»XMIN4,XMAX4#YMIN4-f YMAX4t 3» XYb< 1 ) ♦ 
1 2,PLM(I) ,14*ITAPXG> 

DX = 0 . 

DY=1. 

DO 515 J=2»N02 
K=0 

DX=DX+RCOS 

DY=DY+RSIN 

DU=0. 

DV=1. 

DO 510 I=2»N02 
K=K + 1 

DU=DU+RC0S 

DV=0V+RSIN 

A1H0RZ(K)=J+DU 

A1 VERT (K ) = SpHIM' tt A4 < J» I ) ♦DV 

A2H0RZ (K ) =1 +DX 

510 A2VERT(K)=SPHIM*A4(I»J) *DY 

CALL ODIPLTtO* INGtK»A2H0RZ,A2VERT*XMlN4>XMAX4tYMlN4*YMAX4» 

1 3.XYL(1) *2»PLM(1) *14,ItAPXG) 

CALL DDIPLTtO* ING.K,AlH0RZ,AlV£RT*XMlN4»XMAX4*YMlN4»YMAX4t 
1 3,XYL<1) *2,PLM(1) ,14,ITAPXG) 

515 CONTINUE 

K=0 

DO 520 I-2*N02 
XI = I 

K=K + 1 

A2H0RZ (KJ=Xl+RCOS 
A2VERT<K)=I.+RSIN 
520 CONTINUE 

CALL DDIPLT ( 1 ♦ ING»K* A2HQRZ * A2VERT »XMIN4* XMAX4* YMIN4* YMAX4» 


\ 

V 



1 3, XYL { 1 ) ,2»PLM(1) *13,ITAPXG) 

C HAKE A CONTOUR PLOT OF THE TEMPERATURE 

CALL DDIPLT (0*ING*5*XPR»YPR*XMIN4-,XMAX4,YMIN4,YMAX4*3»XYU(1) * 
1 2»PLT (1) *14* ITAPXG) 

DX=0. 

DY= 1 . 

DO 540 J-2»N02 
K=0 

DX=DX*RCOS 

DYsOY+RSIN 

ou=o. 

DV=1 • 

DO 530 I=2*N02 
K=K + 1 

DU=QU+RCOS 

qv=dv*rsin 

A1H0RH(K)=:J+DU 
A 1 VERT {K)=SPHIT*A3 (J»I ) +DV 
A2H0R2 (K)=I>DX 

530 A2vERT<K)=SpHlT*A3(I»J)+DY 

CALL DDIPLT (0»ING*K»A2HORZ» A2VERT*XMIN4»XMAX4#YMIN4»YHAX4» 

1 3,XYL{1) ,2,PLT(1) ,14, ITAPXG) 

CALL DDIPLT ( 0 ♦ ING*K* A1HORZ » A1 VERT » XMIN4* XMAX4* YM IN4* YMAX4* 

1 3,XYL(1) »2*PLT<1> *14, ITAPXG) 

540 CONTINUE 
K=G 

DO 545 I=2*N02 

XI-I 

K=K + 1 

A2HORZ {K)=XI+RCOS 
A2VERT(K)=I.+RSIN 
545 CONTINUE 

CALL DDIPLT ( 1 * ING»K» A2HORZ» A2VERT *XMIN4*XMAX4» YMIN4* YHAX4* 

1 3, XYL(1)»2«PLTU)*13, ITAPXG) 

C HAKE A CONTOUR PLOT OF THE PRESSURE 

CALL DDlPLT(0,ING,5»XPRtYPR»XMIN4,XMAX4*YMIN4*YMAX4*3»XYU(l) ♦ 
1 2.PLP(1) *14»ITAPXG> 

DX=0. 

DY=1. 

DO 560 J-2*N02 


K = 0 

dx=dx+rcos 
DY=DY*RSIN 
DU = 0. 

DV=1. 

DO 550 1=2, N02 
<=K + 1 

DU=DU+RCOS 

dv=dv+rsin 

A1H0RZ ( K ) = J+DU 

AlVERT(K)=SpHIP*A3(J»I>*A4< j,i>+dv 
A2HORZ (K ) =1 +DX 

550 A2vERT(K)=SpHlP*A3< I*J)*A4(I, J) +DY 

CALL DDIPLT (0* ING,K,A2H0RZ,A2VERT,XMIN4,XMAX4*YMIN4*YMAX4, 

1 3,XYL(1) ,2,PLP(1) ,14,ITAPXG) 

CALL DDIPLT (0, ING,K » A1HORZ , A1 VERT ,XMIN4» XMAX4, YM IN4, YMAX4, 

1 3,XYLU> ,2,PLP(1) ,14,ITAPXG) 

560 CONTINUE 
K=0 

DO 570 1=2, N02 
X 1 = 1 
K=K + 1 

A2HORZ (K)=XI+RCOS 
A2VERT(K>=I.+R5IN 
570 CONTINUE 

CALL DDIPLT (l,lNG,K,A2H0RZ,A2VERT,XMlN4,XMAX4,YMlN4»YMAX4t 
1 3 , XYL ( 1 ) ,2 ,PLP (1) »13 »ItAPxG) 

600 IF(CYY.LT.CY) GO TO 650 

IF(CYY.EQ.CY) WRITE INFORMATION TO BE SAVED ON TAPE1 

TRANSFER PARTICLE POSITIONS FROM TAPE JTG ONTO TAPE1 • 

NS2=0 

610 READ(JTG) xhold,yhold 
WRITE { 1 ) XHOLD*YHOLD 
NS2=NS2+NBSg 

IF{NS2.LT.NBRG) GO TO 610 

WRI-TE(I) Al, A2,A3,A4»CYY,SUMMAS,SUMPX,SUMPY,SUMKE,5UMIE 

rewind 1 

650 IF(IPRNTMG.EQ.O) GO TO 700 
PRINT 652,SCALM 

652 FORMAT (lHl^THE FOLLOWING IS A PLOT OF MASS DENSITY (ALREADY SCALED 


1 WITH DT SQUARED) WHICH HAS BEEN DIVIDED BY *E16.8* AND THEN*) 
PRINT 653 

653 FORMAT (1H *INTEGERI ZED, WITH VALUES BELOW 1 OR ABOVE 9 REPLACED BY 

l blank or star respectively*/) 

DO 670 1=1, N02 
DO 670 J=1»N02 
IM=A4(I,J)/SCALM 
IF(IM.EQ.O) GO TO 660 
IF(IM.GT.9) GO TO 665 
ENCODE (10»655»A4(I»J) ) IM 
655 FORMAT (I 10) 

GO TO 670 
660 A4 { I * J ) -1 OH 

GO TO 670 

665 A4(I»J)=10H * 

670 CONTINUE 

DO 672 1=1, N02 
A4(I»1)=10H 
A4 ( I ,N02) =10H 
672 CONTINUE 

DO 674 J=2*N02M1 
A4 < 1 * J ) =10H I 

A4 (N02 , J ) =1 OH I 

674 CONTINUE 

DO 693 JJ=1 , N02 
J=N02P1-JJ 

693 PRINT 695, <A4<I,J) »I=1,N02> 

695 FORMAT (1H ,64R2) 

PRINT 697 

697 F0RMAT(1H1*IH ) 

700 RETURN 
END 


a,; 



OVERLAY (IFILE»4»Q) 
PROGRAM GASPLOT 


COMMON/ALLCOM/I2A, I TEST *N*CY,CYY»R HO (33*33) 

COMMON/GASCOM/NPLOTG,NPRINTG,NPRNTMG»N02SQ,N02,N02P1,N02MI,DMG* 
1 GM102*AU02*AU04,SavPE(200) fSAVlE(200) *SAVKE(200) tSAVTE(200). 


2 CYMTN.NBRG.NBSG»JTG»jSG»NMKG»PI,MA 5K1*ING<2) »XMlNG,XMAXGf 

3 YmING* YMAXG»XPG( 2) ,YpG, lTAPXG»XYL (3) »PLP 12) »PLT <2) »PLM (2) * 

4 RC0S,RSlN f XPR(5) ,YPR (5) , XMAX4-, XM IN4 , YMAXA, YM IN4, SPHlM* SPKIT, 

5 SpHIP»EMIN*EMAX«ScALM»DT»DT2»DT3*DT4*SUMMAS t SUMPX»SUMPVf 

6 SUMKE*SUHIE*SPHI»PL (2) 


DIMENSION XDATA(200) »LABLEN<8) 

DIMENSION ENDCYY (2) »ENDPE (2) »ENDIE(2) *ENOKE (2) »ENDTE(2) 


INTEGER CY 


CALL PSEUDO 


CYMAX=CY 
NCYY=CY-CYMIN 
DO 20 1=1 »NCYY 
20 XDATA(I)=CYMIN*I 
LABLCYY=10HCYCLES 
L'ABLEN ( 1 ) = 10HENERGY***C 
LABLEN(2)=10HIRCLE-POTE 
LABLEN(3>=10HNTIAL*SQUA 
LABLEN (4 )=10HRE- INTERNA 
LA8LEN (5 ) =10HL ,DIAMQND- 
LABLEN(6)=10HKINETIC»TR 
LABLEN (7) =10HI ANGLE-TOT 
LABLEN (8)=10HAL 
ENDCYY(1)=CYMIN+1 


ENDCYY (2> =CY 
ENDPE(1)=SAVPE(1) 

ENDPE(2)=SAVPE(NCYY) 

ENDIE(1)=SAVIE(1) 

ENDIE<2)=SAVIE(NCYY) 

ENDKEU)=SAVKE(1) 

ENDKE(2)=SAVKE(NCYY) 

END-TE ( 1 ) =SAVTE < 1 > 

ENDTE(2)=SAVTE(NCYY) 

IF (NCYY • LT ■ 2) GO TO 40 

CALL DDIPLT(0,INGt2, ENDCYY, ENDPE*CYMIN 9 CYMAX. EMIN, EMAXt 
1 1 ,LABLCYY, 0, LABLEN, 1 ) 



i :V 

: £ 
■% 






I i ? . 
lr; 



nr 

il 



L'l 



CALL DDIPLT (0» ING» 2» ENOCYY ^ENDI E*CYM IN*CYMAX* EMIN* EM AX t 
1 1,LABLCYY*8*LABLEN*2} 

CALL DDIPLT ( 0 * ING* 2* ENDCY Y * ENDKE*CYM IN *CYMAX * EMIN*EMAX* 
1 l,LABLCYY,8fLABL£N»3) 

CALL DDlPLTlO*ING*2fENDCYY,ENDTE*CYMlN*CYMAX*EMIN*EMAX* 


1 1,LABLCYY,8,LABLEN,4) 

40 CALL DDlRLT(0*ING»NCYY*XDATAfSAVPEtCYMIN#CYMAX*EMIN*EMAX* 
1 1 tLABLCYY ♦ 8»LABLENt 0 ) 

CALL DDIPLT ( 0 » ING*NCYY * XDATA* SAV IE» CYM INtCYHAX* EM IN*EMAX* 
1 1*LABLCYY*8*LABLEN»0) 

CALL DDIPLT (0»ING*NCYY»XDATA,5AVKEfCYMIN»CYMAX»EMIN,EMAX* 
1 1 ,LABLCYY#8tLABLEN» 0) 

CALL DDIPLT ( 1 » ING*NCYY» XDATA* SAVTE* CYM IN* CYM AX* EM IN* EM AX* 


1 1,LABLCYY,8*LABLEN*0) 

RETURN 

END 



APPENDIX E 


Computer Plots Produced by the Jeans Instability Simulator of 
Appendix D. 

These plots are from two full length runs with identical initial 
conditions except for the ratio of total internal to total potential 
energy. The mesh size is 32 x 32. Later runs will be made with meshes 
of 64 x 64 and possibly 128 x 128. 

The last plot of each run demonstrates that despite rapid changes 
in total potential, internal, and kinetic energies, the code conserves 
total energy almost exactly. 
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APPENDIX F 


Listing of the Two-Dimensional Polar Coordinate Particle-in-Cell 
Simulator of a Rotating Self-Gravitating Compressible Gas. 

Comment cards necessary to make this listing self explanatory will 

be added later. 


Overlay No . 

Program/Subroutine Name 

Page No. 

Co ,o) 

PROGRAM POLR 

F-2 

U ,o) 

PROGRAM GETPHI 

F-3 

C2 ,0) 

PROGRAM INIGAS 

P-5 

(3,00 

PROGRAM ADVGAS 

F-16 


SUBROUTINE NEGIE 

F-35 

C4 ,0} 

PROGRAM GASPLT 

F-41 






OVERLAY (IFILE,0,0> 

PROGRAM POLR (OUTPUT, TAPE] , rAPEB * JAPE? , TAPE 1 0 » TAPE 1 1 ) 
COMMON/ALLCOM/I2A, I TEST , N , N02 , CY , CYY * RHO (32,32) 

COMMON /GAS COM /NTH (15) , DTH (15) , AREA (15) *NPL0t6,NPRNTG,NRING,N04, 

1 N04M1,N04M2,RCC,GMC,GMCP1,GMCP2,DMG*DT , DT2 » DT3, DT4, NRTCEL', 

2 RMAX * CENTER, SUMPTH.SUMM AS, SUM IE ,SUMKE, RESOLD* PEFOLO»XMG* 

3 SAVPTH(200) »5avPE(200) ,SAVIE(200) ,SAVK£(200) ,SAVTE<200) * 

4 MASK1 ,NBSG,NBRG,JTG, jSG, IMG(2) , XMING»XMAXG, YM ING, YMAXG, XPG ( 2 ) 

5 YPG»ITAPXG,GM1,NMKG,XYL(3) , XPR (5 ) , YPR ( 5 ) ,XMIN4,XMAX4,YMIN4, 

6 YMAX4,PLM (2) »RCOS»RSiN»SPHIM,PLT (2) »SpHIT,PLP (2) »5PHlP»TrtOPI » 

7 PTHMlN,PTHMAX»EMlNtEMAX*CYMIN,PI*PL(2 ) ,SPHI»MTESt * I TAPE 


INTEGER CY,CYY - 
IFILE=5LIFILE 

GPF ILE=6LGPF ILE 
AGFILE=6LAGFILE 

mtest=i 

itegt=i 

CALL OVERLAY(IFIlE, 2,0,6 HRECALL) 

IF ( I T APE »EQ d 1 ) CALL OVERLAY (GPFILE, 1 * 0 » 6HRECALL-) 
IF ( ITAPE.EO® 1 ) CALL OVERLAY ( IF ILE » 2 , 0 * 6HRECALL’) 

5 PRINT 10, CYY 
10 FORMAT (1H0»CYCLE = *18) 

CALL SECOND (Tl ) 

CALL OVERLAY (GPFILE, 1 , 0 , 6HRECALL ) 

CALL SECOND ( T2 ) 

T3=T2-T 1 
PRINT 15, T3 

15 FORMAT ( 1 H *FIELD TIME = *E16.8) 

CALL OVERLAY (AGF ILE, 3,0* 6HRECALL) 

CALL SECOND (T4) 

TS=T4-T2 
PRINT 35sT5 

35 FORMAT ( 1 H ^ADVANCE TIME = *E16.8) 

IF (CYY.GE.CY) GO TO 40 


CYY=CYY ♦ l 
GO TO 5 

40 CALL OVERLAY ( IFILE*i4* 0 * 6HRECALL ) 
STOP 
END 




OVERLAY (GPFILE. 1.0) 

PROGRAM GETPHI 

C0MM0N/ALLC0M/I2A» iTESTf N» N02»CYf CYY*RHOA(32*32) 
COMMON Z(1025)*Y(1025) 

DIMENSION G(33.33> .RH0<64. 32) 

I2B=I2A-1 

N21=N02*1 

IFUTEST.EQ.O) GO TO 10 

ITEST=0 

RNI = 1 •/ (N*N) 

DO 1 J=1 f N2 1 
DO 1 I=ltN21 

IF(I.EQ.I.AND.J.EQ.l) GO TO 1 

G( If J)=RNI/SQRT( < I-1.)*(I-1.) ♦ < J-l . )*( J-l .) ) 

1 CONTINUE 

G( 1 » 1 ) =G ( 1 «2) 

CALL GETSET(2f I2B) 

DO 2 J=1 . N2 1 
DO 3 1=1. N21 

3 Z(I)-G(I»J) 

CALL FTRANS<2» I2B) 

DO 4 1=1 t N21 

4 G { I f J ) =Y 1 1 ) 

2 CONTINUE 

DO 5 1=1 tN21 
DO 6 J=1 . N21 

6 Z(J)=G(I*J> 

CALL FTRANS(2*I2B) 

DO 7 J=1.N21 

7 G< I . J) =Y ( J) 

5 CONTINUE 
WRITE ( 1 ) G 

REWIND 1 
10 READ { 1 ) G 
REWIND 1 

CALL GETSET(3.I2A) 

DO 11 J=1 . N02 
DO 8 1=1 . N02 
Z ( I ) =RHOA ( I » J > 

8 Z (N02* I ) =0 • 


CALL FTRANS (3, I2A) 

DO 9 1 = 1, N 
9 RHO ( I , J ) =Y ( I ) 

11 CONTINUE 

DO 12 1=1, N 
DO 13 J=1,N02 
Z ( J) =RHO ( I » J ) 

13 Z ( J+N02 ) =0 • 

CALL GETSEH3,I2A) 

CALL FTRANS(3,I2A) 

IFU.GT.N21) GO TO 14 
DO 15 J=2, N02 
Z < J) -Y ( J) tt G{ I * J) 

15 Z(J+N02>=Y(J*N02)*G<I*J> 
Z(1)*Y(1)*GIM) 
Z(N21)=Y<N21>*G(I,N21) 

GO TO 16 

14 DO 17 J=2, N02 
Z<J)=Y<J)*G<I-N02,J) 

17 Z ( J+N02) =Y ( J+N02) *G ( I-N02* J) 
Z(1)”Y(1)*G( I -NO 2, 1 ) 

Z<N21)=Y(N21)*G<I-N02,N21) 

16 CONTINUE 

CALL GETSET(4,I2A) 

CALL FTRANS (4, I2A) 

DO IB J=1,N02 

18 RHO ( I » J) =Y ( J) 

12 CONTINUE 

DO 19 J=1,N02 
DO 20 1=1, N 

20 Z ( I ) =RHO ( I , J) 

CALL FTRANS (4, I2A) 

DO 21 1=1, N02 

21 RH0AU,J)=Y(I> 

19 CONTINUE 

return 

END 










s-*-~ «*»'«' i'i- n 


OVERLAY (IFILE»2»0) 

PROGRAM INIGAS 

COMMON/ ALL COM/I 2 A* I TESj« N* NQ2 * CY , CYY * RHO ( 32* 32) 

COMMON/GASCOM/NTH ( 15) ♦ DTH ( 1 5 } ,AREA(15) » NPLOTG , NPRNTG * NR I MG* N04 ♦ 

_ _ .. . . _ * — n r-% . . .r* r-k— 1. 1 1 H T P C I . 


SET 


N04M1 »N04M2*RCC*GMC*GMCPi*GMCP2,DMG*DT * DT2 * DT3 * DT4* NRTCED* 

RMAX*CENTER*SUMPTH,SUMMAS,SUMIE»SUMKE»PESOLD*PEFOLO*XMG, 

5avPTH(200) ,SAVPE<200) «SAvIE(200) *SAVKE (200 ) *SAVTE(200)» 

MASK1 tNBSG»NBRG* JTG* JSG* ING ( 2) ♦ XM ING* XMAXG* YM ING, YMAXG, XPG ( 2 > ♦ 
YPG* ITAPXG,GM1 *NMKG, XYL (3) * XPR ( 5 ) »YPR (5) »XMIN4»XMAX4*YMIN4, 
YMAX4*PLM (2) ,RCOS,RSlN*SpHlM*PLT (H) »SPHIT »PLP (2) .SPHIP.TWOPI* 
PTHMIN»PTHMAX*EMIN*EMAX*CYMIN»PI*PL(2) .SPHIfMTESr » IT APE 
DIMENSION RADIUSdOOO) , THET A ( 1000 ) •ROOT (15,64). TOOT (15,64). 

1 IE (15*64) *MASS< 15,64) , PRES (15,64) tXNWAV ( 15*64) *YNWAV ( 15.64) « 

2 R2MDAV (15*64) ,R2NWAV ( 15*64) *GR(15> ,KAPPA2(15) 

INTEGER CY.CYY 

REAL IE*MINIE»KAPPA*KAPPA2. MASS, NUMBER 
INITIAL VALUES 

SET I2A WHERE (2#*I2A)/2 IS THE SIZE OF THE ACTIVE MESH 
I2A-6 

SET NUMBER OF CYCLES 
CY=60 

SET TOTAL NUMBER OF GAS PARTICLES 

NBRG=1000~ , _ 

GAS PARTICLES PER READ OF GAS PARTICLE FILE 


C 

C 


i c 


set numb: 

NBSG=IO r 
SET INIU 
RIG“ 1 0 • 

SET TOTAL 
GMC=3.55E6 

SET FRACTION OF GALAXY MASS WHICH IS GAS 
PERC= • 05 
SET RADIUS OF 
RCC=6. 

SET INITIAL MAXIMUM 
WIDTH. 

ARATlO= • 3 ^ _ 

SET NUMBER OF POINTS PER GAS PLOT 

npg=nbrg 

SET PERIOD OF GAS LONG PRINTING 

SET PERIOD OF PRINTING OF GAS ANGULARLY AVERAGED (RING) QUANTITIES 


„ RADIUS OF GAS 
MASS OF GALAXY 

" GALAXY MASS WH 

fixed sphere of stars 

ratio of time scaled angular 


velocity 


to ANGULAR CELL 


| MR ING=30 

| C set period OF gas PLOTTING 

| c set°ratio of specific heat at constant volume to specific heat at 

| c constant pressure, (gamma must be greater than or equal to one. 

I C GAMMA IS EQUAL TO 2.0 FOR NORMAL SIMULATION OF MONATOMIC GAS IN 

| C TWO DIMENSIONS. GAMMA MAY BE SET EQUAL TO 1.0 FOR A SIMULATION 

I c WITHOUT PRESSURE TERMS.) 

:| c SET H RATI0 OF THE VELOCITY DISPERSION OF THE GAS (SORT OF SPECIFIC 

I C INTERNAL ENERGY) TO THAT REQUIRED FOR LOCAL STABILITY 

QGAS=1.0 

! C SET MINIMUM INITIAL SPECIFIC INTERNAL ENERGY 

I MINIE=1. 

I c SET ITAPE=1 TO START RUN. SET ITAPE=2 TO CONTINUE RUN WITH PICK UP TAPE, 

f ITAPE=1 

| C SET CONSTANTS 
4 PI=3 • 1415926536 

TWOP 1=2 .*P I 

MASK 1=07777777777000000 0000 

| N=2**I2A 

N02=N/2 
f N04=N/4 

I N04M1=N04-1 

j N04M2=N04-2 

% NRTCEL=N*N04M1 

f JTG=9 

I JSG=10 

! GM1=GAMMA-1 • 

RIG2=RIG**2 
NMKG=NPG-NBSG 
j XMG=GMC*PERC/NBRG 

f XMG=XMG.AND.HASK1 

I GMC=GMC*(1.-PERC) 

4 CENTER=N04v,5 

I RMAX=N04“2 .0001 

| NTH ( 1 ) =4 

|: NTH (2) =8 

NTH ( 3 ) =16 
n NTH (4) =16 



j 


1 


■% 


li 



NTH(5> =32 
NTH (6) =32 
NTH ( 7 ) =32 
NTH ( 8 > =32 
NTH ( 9 ) =64 

NTH (10) =64 
NTH ( 1 1 ) =64 
NTH (12) -64 

NTH { 13 ) =64 
NTH ( 14 ) =64 
NTH (15) =64 

IF (N04.EQ ■ 16) GO TO 20 
NTH (16) =64 
DO 10 I = 17» 31 
NTH ( I >=128 
10 CONTINUE 
20 DO 30 1=1 »N04M1 

DTH ( I ) =T WOP I /NTH ( I ) 
AREA(I>=DTH<I)*(I-.5> 

30 CONTINUE 

C SET PLOTTING SPECIFICATIONS 
ING ( 1 ) =1 0H2D-GAS 
ING(2)=10H 
XMING=-10 
XHAXG=40 
YMING=-10 
YMAXG=40 

XPG<l)=10HXtCYCLE= 

YPG= 1 OH Y 
IT APXG=6LT APE23 
XYL ( 1 ) =1 OHX-Y PLANEt 
XYL (2 ) =1 OH CYCLE = 
PU1>=10HP0 tENTIAL* 

PLP ( 1 ) = 1 OH PRESSURE* 
PLTll)=10H TEMP* 

PLM ( 1 ) =1 OH DENSITY* 
ANG=PI/3. 

RCOS=COS ( ANG) 
RSIN=SIN(ANG) 

IC=N02+1 


XPR ( 1 ) =1 • +N02#RC0S 
XPR(2>=1. 

XPR(3)=1.*N02 
XPR C ^ ) =IC*N02*RC0S 
XPR { 5 ) =XPR ( 1 ) 
YPR<1)=1.*N02*RSIN 
YPR(2)=1. 

YPR (3) =1 • 

YPR(4)=1.+N02*R5IN 

YPR(5)=YPR(1) 

XMAX4=IC+N02*RCOS 
XMIN4=-1 . 

YM IN4=-1 • 

YHAX4=XMAX4 

HT = ION02«RC0S-N02*RSlN 
END OF INITIAL DATA 

IF ( I T APE .EQ.2) GO TO 350 
CYY= 1 

IF{MTEST.EQ.O> GO TO 160 
MTE5t=0 

DO 105 1=1 »N04M1 
NNTH=NTH(I) 

DO 105 J=1*NNTH 
ROOT ( 1 1 J ) =0 • 

MASS ( I » J) =0 • 

XNWAV( I*J)=0. 

YNWAV { I * J> =0 • 

R2MOAV { I ♦ J) =0* 
R2NWAV(I»J)=0. 

105 CONTINUE 

DO 107 I=1«N02 
DO 107 J=1»N02 
RHO(I f J)=0. 

107 CONTINUE 

X=URAN (7654321.) 

NSa=o 

110 DO 130 IS=1,NBSG 
120 X=2. tt (UR AN (0 • ) -.5) 
Y=2.*(URAN(0.}-.5) 

Z7=2 (URAN ( 0 • ) -.5) 


R2=X tt X+Y*Y+Z7*Z7 
IF(R2.GT.l.) GO TO 120 

XC=RIG*X 

YC=RIG*Y 

IFtXC.EQ.Q.) XC-, 00001 

R2=XC*XC*YC*YC 

R=SQR T (R2) 

TH=ATAN2(YC,XC) 

IF (TH.LT »0 • ) TH=TH*TWCPI 
IF(TH.GE.TWOPI> TH=TK-TWOPl 
IR=R+ 1 

JT=TH/DTH( IR) 4-1 

ix=center*xc*.s 

JY=CENTER*YC«.5 
XNWAV(IR*JT)=XNWAV(IR* JT) + xc 
YNWAV(IRtJT) =YNWAV ( IR*JT) +YC 

R2MDAV ( IR» JT) =R2M0AV ( IR» JT> +R2 
R2NWAv(IR»JT)=R2NWAy(IR*jT) +R2 
MASStlR, JT)=MASSaR*JT> +XMG 
RHO(IX,JY)=RHO(IX,JY> *.5»XMG 
RADIOS ( IS) =R 
THETA (lS)=TH 
130 CONTINUE 

WRITE (9) RADIUS, THETA 
NS2=NS2+NBSg 

IF(NS2,LT.NBR6) 60 TO 110 
REWIND 9 

DO 1A0 I=1»N0AM2 
NNTH=NTH(I ) 

DO 1 AO J=1,NNTH 

IF(MASS II, J) .EQ.O.) GO TO 1 AO 

NUMBER=MASS{ I, J)/XMG 
XNWAV ( I * J ) =XNW AV (I * J > /NUMBER 
YNWAV* I t J) =YNWAV ( I , J) /NUMBER 
R2MDAy ( I » J) =R2MDAv 1 1 * J) /NUMBER 
R2NWAV < I , J ) =R2NWAV ( I * J) /NUMBER 
1A0 CONTINUE 

WRITE (8) ROOT »HASS, XNWAV, YNWAV»R2MDAV,R2NWAV 

REWIND 8 
RETURN 


160 READ (8) RD0T*MASS,XNWAV*YNWAV*R2MDAV»R2NWAV 
REWIND 8 

DO 170 I-UN04M2 
GR ( I ) =0 * 

NMtH^NTH { I ) 

DO 165 J=1»NNTH 
R= I” * 5 

TH= ( J-.5) S DTH ( I ) 

XCELL=R*COS(TH) 

YCELL=R*SIN(TH> 

X=CENTER+XCELL 

Y=CENtER+YCELL 

lx=x 

I Y=Y 

xx=ix-x 

YY=I Y-Y 
RAD3=R*R*R 

IF(R.LT.RCC) RAD3=RCC*RCC*RCC 

GX=(YY*1.)*< (XX+l a )*(RHO(lX+l*lY)-RHO(IX-l*lY) ) “XX* (RHO (IX*2»IY)-R 
1 HO ( I X * I Y ) ) ) - YY* ( ( XX*1 • ) * (RHO ( IX+1 * IY+ 1 )-RHO (IX-1*IY*1) ) -XX* (RHO (IX 
2+2»IY*l) -RHO ( IX ♦ I Y* 1 ) ) ) -6MC*XCELL/R AD3 
GY=(XX + 1 . >*< ( Y Y + 1 • ) * ( RHO ( I X * I Y + 1 ) -RHO ( IX ♦ I Y- 1 ) ) -YY* ( RHO < IX * I Y*2 ) -R 
1HO { IX* IY ) ) >-XX*( (YY«-1.)*(RH0(IX*1*IY-*1>-RH0(IX*1*IY-1) ) - YY* ( RHO ( I X 

2+1iIY+2)-RH0(IX*1*IY> ) ) -GMC* YCELL/RAD3 
GR ( I ) =GR (!)•*•( XCELL*GX + YC£LL*GY ) /R 

165 CONTINUE 

GR ( I ) =GR ( I ) /NNTH 
170 CONTINUE 

N04M3=N04~3 
DO 175 I=2*N04M3 
R= I- . 5 

KAPPA2(I)=(GR(I+1)-GR(I-1) )/2.+3.*GR(I)/R 
175 CONTINUE 

KAPPA2U )=(GR(2)-0.)/2.+3.*GR(l)/.5 

KAPPA2 (N04M2) = (GR (N04M2) -GR (N04M3) ) /I • ♦3«*GR ( N04M2) / ( N04M2- *5) 
KAPPA 2 (N04M 1 ) = 0 • 

DO 185 1=1 »N04M1 
KAPPA2 ( I ) =-K APPA2 ( I ) 

I F (K APPA2 ( I ) .LT • 0 » ) KAPPA2(I)=0. 

KAPPAsSQRT (KAPPA2 ( I > ) 

fl NNTH=NTH(I) 


DO 185 J=1*NNTH 
IF(MASSdtJ) .EQ.O.) GO TO 182 
IF (KAPPA. EQ.O. ) GO TO 178 

IE(It J)= (3.36*QGAS*MASS< I, J ) / ( AREA d ) *KAPPA ) )**2+MINIE 
GO TO 180 
178 IE (I ♦ J) =MINIE 

180 PRES( I , J)=GM1*IE( It J»*MASS( I , J) /AREA ( I ) 

GO TO 185 
182 IE ( I « J ) -0 • 

PREG ( 1 1 J ) =0 » 

185 CONTINUE 
AMAX=0. 

ISAV=1 

DO 240 1=1 tNQ4M2 
IPTESt=0 

IFd.EQ.l) IPTEST=1 
IFU.EG.2) IPTEST=1 
IFd.EQ.4) IPTEST=1 
IF(I.EQ.8) IPTESt=1 
IF (I.EQ.16) IPTEST=1 
IMTESt=0 

IFU.EQ.2) IMTEST = 1 
IF ( I .EQ.3) IMTEST=1 
IF(I.EQ.S) IMTEST=1 
IF(I.EQ.9) IMTEST-1 
IFd.EQ.17) IMTESt=1 
NNTH=NTH(I> 

DO 240 J=1 tNNTH 

IF(MASS(It J) .EQ.O.) GO TO 230 

XNEW=XNWAV(ItJ) 

YNEtf-TNWAVdt J) 

RNE*=SQRT (XNEW*XN£W+YNEW»YnEW) 

X=CENTER+XNEW 

y=center*ynew 

I X=X 
IY=Y 
XX=IX-X 
YY=IY-Y 

RAD3=RNEW*RNEW*RNEW 
IF(RNEW.LT.RCC) RAD3=RCC*RCC*RCC 

GX=<YY*1.)*{ <XX«-l.)*(RHOdX + ltIY)-RHOdX-lt IY) > -XX* (RHO <IX*2t IY) -R 


1 HO ( IX, I Y) ) )-YY*( (XX* I.) *<RHO( IX+1* IY*1)-RH0( IX-1» IY + 1 ) >-XX* (RHOtlX 

2*2*IY+1)-RH0(IX*IY+1) ) ) -GMC*XMEW/RAD3 
GY=(XX*1.>*( ( YY * X • ) * ( RHO ( I X ♦ I Y * 1 ) -RHO (IX*IY-1) ) -YY* (RHO ( IX * I Y*2 ) -R 
1 HO ( IX » I Y > ) )-XX*< (YY*l.)*(RHO(IX*l* I Y* 1 ) -RHO ( I X* 1 ♦ IY-1 ) >-YY*(RHO(lX 
2*1»IY+2)"RH0(IX*1* IY) ) ) -GMC*YNEW/RAD3 
GRR=(XNE**GX+YNEW*GY)/RNFW 

IF(IPTESt.EO.I) go TO 200 

PIP=(PREG{I,J) + PRESU + 1,J> )/2. 

IF (MASS ( 1 + 1 * J) .EQ.O • ) PIP = 0. 

GO TO 210 

200 PIPI=(PRES(I,J)*PRES(I*l t 2*J-l) )/2. 

IF (MASS (1+1.2* J-l) .EQ.O.) PIP 1 = 0. 

PIP2=(PRES(l,j) +PRES ( I + 1 . 2*J ) )/2. 

IF (MASS(I+1.2*J) .EQ.O.) PIP2=0. 

PlP”(PIPi+PIP2)/2. 

210 IF(I.NE.l) GO TO 215 
JP2- J+2 

IF ( JP2.GT .A) JP2=JP2~4 
PIM=(PRES(1, J)+PRES{1, JP2> )/2. 

IF(MASS(1,JP2) .EQ.O.) PIM=0. 

GO TO 220 
215 JJ=J 

IF(IMTEST.EQ.l) JJj= ( J+ 1 ) /2 
PIM=(PRES(I, J) *PRES(I-l,jJ) )/2. 

IF(MASS(I-1,JJ) .EQ.O.) PIM=0. 

220 TD0T2=- ( (PIH-PIP) /MASS ( I , J) + GRR) /RNEW 
IF (T DO T2.LT .0.) TDOT2=0. 

TDOT(I»J)=-SORT(TDOT2) 

ASAV=-TDOT ( I * J ) /DTH ( I ) 

IF(ASAV.LT. AMAX) GO TO 240 

amax=asav 

ISAV=I 
GO TO 240 
230 TOOT ( I * J) =0 . 

240 CONTINUE 

NNTH=NTH(N04M1) 

DO 245 J=1*NNTH 
TOOT (N04M1 * J) =0 • 

245 CONTINUE 

dt=aratio/amax 
PRINT 250. DT 


250 format(1ho*time step dt = *Ei6.e> 

DDDG=NTH ( ISAV) /ARATIO 
PRINT 253* ISAV*DDDG 

253 FORMAT(1HO*AT R«-.S = *13* THERE ARE *F1G.3* TIME STEPS PER ROTAT 
I ION* ) 

DT2=DT*DT 

DT3=DT2*DT 

DT4=DT2*DT2 

DMG=XMG*DT2 
DMG=OMG. AND. MASK 1 

PRINT 255 * DMG 

255 FORMAT(IHO*TIME ScALEO PARTICAL MASS = *E16.8) 

GMC=GMC*DT2 

GMCP1=3.*GMC/(2.*RCC) 

GMCP2=GMC/ (2 . *RCC**3 ) 

DO 260 1=1 *N04M2 

NNTH=NTH(I) 

DO 260 J=1*NNTH 
IF (MASS ( I * J) . EQ.O. ) GO TO 260 
TDOT(I*J)=TDOT (I*J)*DT 
IE(I,J)=IE(I*J)*DT2 
MASS ( I , J)=MASS ( I * J) *DT2 
260 CONTINUE 

ANGVG=TWOPI/(DDDG*DT> 

EMAX=.2*XMG*NBRG*RIG2*ANGVG*ANGVG 

EMIN=-10.*EMAX 

PTHMAX=0. 

PTHMlN=~.5*XMG*NBRG*RIG2*ANGVG 
PMAX = 2.*(2.*RH0(N04*N04)*DT2 + GMCP'1> 

PMAXM=10.*DMG*NBRG/(PI*RIG2> 

PMAXT=I00.*(IE(6*1) ♦IE(6*2) ♦ IE (6.3) ♦ IE <6. 4) )/4. 

P MAXP= • 2*PMAXM*PMAXT 
SPHI =HT /PMAX 

ENCODE ( 1 0 *275* PL < 2 ) > SPHI 
275 FORMAT (F10.3) 

SPHIM=HT /PMAXM 

ENCODE ( 10*275 *PLM (2) ) SPHlM 

sphit=ht/pmaxt 

ENCODE ( 10*H75*PLT (2) ) SPhIT 
SPHIP=HT/PMAXP 


ENCODE (10*280 »PLP (2) ) SPHlP 
280 FORMAT(FIO.O) 

summas=o. 

sumpth=o. 

SUMKE=0. 

sumie=o. 

DO 300 1 = 1 *N04M2 
NNTH=NTH< I) 

DO 300 J=1#NNTH 

IF (MASS { I , J ) ■ EQ • 0 . ) GO TO 300 

R2MI0=R2MDAV<I* J) 

SUMMAS=SUMMAS+MASS { I * J) 

SUMPTH=5UMPtH*MAS 5 ( I » J > e R2M ID^TDOT < I * J) 

SUMKE=SUMKE* .5*MASS(I, J)*R2MID*TD0T( I*J)*TD0T (I* J) 
5UMIE=SUMIE*MAS5(I,J)«IE(T*J) 

300 CONTINUE 
GO TO 4-00 

c statements 350 to 400 enable run to be continued with pick up tape 

350 NS2=0 

355 READ ( 1 1 ) RADIUS, THETA 
WRITE (9) RADIUS, THETA 
NS2=NS2*NBSg 

IF(NS2.LT.NBRG) GO TO 355 
REWIND 9 

READ (11) RDOT*TDOTt IE,MASS,XNWAV*YNWAV»R2MDAV»R2NWAV»5UMMAS, 

1 SUMPTH,SUMIE»SUMKE,PES0LD»PEF0LDtCYY*DT*DT2*DT3*DT4*DMG,GMC, 

2 GMCP1*GMCP2«SPHIM*SPHIT*SPHIP.,PLM<2) *PLT (2) ,PLP (2) *EMAX» EMIN* 

3 PTHMAX*PTHMIN*SPHI,PL (2) 

CY=CY+CYY 

CYY=CYY+1 
400 CYMIN=CYY-1 

CALL EVICT (6LTAPE11) 

DO 405 1=1 *N02 
DO 405 J= 1 , N02 
RHOCI, J)=0. 

405 CONTINUE 
NS2=0 

410 READ (9) RADIUS, THETA 
DO 420 IS=1,NBSG 


IF (RADIUS ( IS) .GT.RMAX) GO TO 420 
IX--RADIUS( IS)*C0S(THETA(IS) ) +CENTER+.5 


JY=RAOIUS ( IS) »S IN (THETA ( I S > ) +CENTER+.S 
RHO(lX»JY)=RHO(IX»JY) +*5*DMG 
420 CONTINUE 


NS2=NS2+NBSG 

IF(NS2.LT.N8RG) GO TO 410 
REWIND 9 


WRITE (8) ROOT »TDOT* IE, MASS, XNWAV,YNWAV*R2MDAV»R2NWAV 
REWIND 8 

RETURN 


END 


OVERLAY (AGFILE*3*0) 

PROGRAM ADVGA5 

COMMON/ ALLCOM/ 1 2A » ItESt *N»N02«CY « CY Y * RHO (32*32) 

COMMON /GASCOM /NTH ( 15)*DTH(15)*AREA(15) ♦ NPLOT G * NPRNtG * NR I NG * NO** ♦ 

1 NOAM 1 ,NOAM2,RCC,GMC,GMCPl * GMCP2 * DMG* DT » DT2 * DT3 * OTA * NRT CEL * 

2 RMAX * CENTER *SUMPTH, SUMMAb, SUM IE * SUMK E * PESOLD* PEFOLD* XMG* 

3 SavPTH(200) *SaVPE( 200) *SAVlE(200 ) »SaVKE(200> *SAVTE(200) » 

A MASK1 *NBSg*NBRG* JTG* jSG* I NG ( 2 ) * XM ING* XMAXG * YM ING* YMAXG* XPG ( 2) ♦ 

C YPG.ITAPXG*GM1* NMKG ♦ XYl(3) *XPR(5)*YPR(5) *XMI NA »XMAXA*YMIN4, 

6 YMAXA*PLM (2) ,RCOG,RS IN,5PHIM tPLT < 2) *SPHIT ♦ PLP (2) * SPH I P ♦ T WOP I ♦ 

7 PTHMIN*PTHMAX*EM IN*EMAX*CYMI\1,P1»PL(2>,SPH1 *MTES T * I TAPE 
COMMON /NEGCOM/ IE ( 1 5 » 6 A ) ♦ M ASS ( 1 5 * 6A ) * XNW AV ( 1 5 * 6 A ) 

DIMENSION RADIUS (1000) ,THETA(1O0O) *RDOT (15*64) * 


C IF 


1 

2 

3 

A 

5 


C 

C 


C 

c 


TDOT ( 15*6A) *YNWAV (15*64) »R2MDAV (15*64), 

R2NWAV< 15*64) ,PRES( 15*64) ,RDDOT ( 15 ♦ 64 ) , PRAD <15 , 64 ) * 

PTHETA ( 1 5* 6A ) fETOTAL ( 1 5 * 6 A ) * XPLOT ( 1000) * YPLOT ( 100 0 ) * GR ( 15) ♦ 

RNEWRU5) * DENSR (15) *PRESR(15) ,MASSR(15) , IERU5) ♦ KAPPA (15) » 

XPLOT1 (32) ♦ YPLOTl (32) *XPL0T2 (32) *YPL0T2 (32) 

EQUIVALENCE (RADIUS { 1 ) , XPLOT (1 > ) * (THETA (1 ) ,YPLOT U ) ) 

INTEGER CY * CYY ♦ Q 

REAL IE* INE* IER*KAPPA*KAPPA2*KE*MASS, MAS, MASSR, NUMBER __ r n 

gas plotting is to be done this cycle set iplotg=i* otherwise set iplotg-o. 

IPLOTG=0 

IP (CYY-CYY/NPLOTG*NPLOTG,EQ.O .OR.CYY.EQ.l ) IPLOTG=l 
ENCODE (10*15* XPG(2) ) CYY 
ENCODE(10*15*XYL(3) ) CYY 

15 FORMAT(IIO) 

CALL PSEUDO 

IF GAS LONG PRINTING IS TO BE DONE THIS CYCLE SET IPRNTG=1* OTHERWISE. 

SET I PRNTG=0 
IPRNTG=0 

IF (CYY-CYY/NPRNTG*NPRNTG.EQ.O .OR .CYY.EQ. 1 ) IPRNTG=1 
IF GAS RING AVERAGES ARE TO BE PRINTED THIS CYCLE SET IRING=1 *OTHER*lSE 
IRING=0. 

IRING=0 

IF (CYY-CYY/NR ING*NR ING. EQ.O. OR. CYY.EQ. 1 ) IRING=1 
READ (8 ) RDOT ♦ TDOT * IE* MASS* XNW AV » YNWAV ♦ R2MDAV * R2NWAV 

REWIND 8 


SET 


IF(IPLOTG.NE.l) GO TO 40 

C MAKE A CONTOUR PLOT OF THE GRAVITATIONAL POTENTIAL 

CALL DDIPLT (0 * ING*5*XPR* YPR*XM IN4« XMAX4* YM IN4* YMAX4*3* XYL't 1 ) • 
1 2,PL(1) »14*ITAPXG) 

DX=0. 

OV= 1 • 

DO 30 J=2iN02 
K=0 

DX=DX+RCOS 

dy=dy+rsin 

DU=0. 

DV = 1. 

DO 25 I=2f N02 
K=K ♦ 1 

DU=DU*RCOS 

dv=dv+rsin 

XC=I-CENTER 

YC=J-CENTER 

R2=XC»XC*YC»YC 

R=SQRj(R2) 

EGC=GMCP1-GMCP2*R2 
IF(R.GT.RCC) EGC=GMC/R 

XPLOT1 (K) =J+DU 

YPLOT1 (K)=SPHI* (2.*RHO t J,I) + EGC) +DV 
XPLOT2 (K ) =I+DX 

25 YPLOT2 (K) =SPHI» (2.*RH0 ( I » J> +EGC) +DY 

CALL ODlPLT(0t INGiK*XPL0T2*YPL0T2»XMIN4»XMAX4*YMIN4»YMAX4» 

1 3,XYLU> *2*PL(1) *14.ITAPXG) 

CALL DDIPLT (0, ING»K»XPLOTl *YPL0T1*XMIN4»XMAX4*YMIN4*YMAX4» 

1 3*XYL(1) *2»PL ( 1 ) *14* ITAPXG) 

30 CONTINUE 
K=0 

DO 35 I=2» N02 
XI = I 
K=K + 1 

XPLOT2(K)=XI*RCOS 









35 


YPL0T2(K)=1.*RSIN 
CONTINUE 

CALL DDIPLT ( 1 * ING»K *XPLOT2* YPL0T2*XMIN4*XMAX4* YMIN4*YMAX4* 

1 3»XYL d ) *2»PL d ) *13*ITAPXG) 

40 IF (IRING.NE.l) GO TO 120 
DO 50 1=1 *N04M2 
GRR=0 . 

NNTH=NTH(I > 

DO 45 J=1 * NNTH 
R=I-.5 

TH={J-.5>*DTH(I) 

XCELL=R*COS(TH) 

YCELL=R*SIN(TH> 

x=center*xcell 

Y=CENtER+YCELL 

ix=x 

I Y=Y 

xx=ix-x 

YY=I Y-Y 
RAD3=R*R*R 

( (Xx2l!l« <RHO?lS“, I Y ) — SHO ( IX-1 . IY> I -XX* (RHO( IX-2. IY) -R 

1HO ( IX* IY ) ) ) -YY* ( (XX* 1 • > ° (RHO ( IX* 1 * IY+ 1 ) -RHO ( IX-1 1 IY* 1 ) )-XX & (RHO(lX 

2 SY=u;UT?*UYY’lIU(RHOmtmn-R5oUX.IY-l>)-YY*IRHO(IX.IY.2)-R 
lHOdXdYdl-XX^dYY + l.l^tRHOdX + DIY + l) - RHO dX + l*IY~l) )-YY*(RH0dX 
2* 1 » I Y+2 ) -RHO (lX+lfIY))) -GMC ft YCELL/RAD3 
GRR = GRR + (XCELL*GX*YCELL*GY ) /R 
CONTINUE 

GR ( I ) =GRR/ <NNTH*DT2> 

CONTINUE 
PRINT 52 

FORMAT <1H0*ANGULARLY AVERAGED QUANTITIES for THE END OF THE LASt C 
IySlE WHICH HAVE BEEN COMPUTED AT CELL GEOMETRIC CENTERS*) 

FORMAT UHO* I R GR KAPPA(-GR)*) 

DO 60 1=1 *N04M2 
R = I - • 5 

IFd.EQ.DG0 TO 54 
IF ( I • EQ .N04HH) GO TO 56 
KAPPA2= (6R(I*1)-GR(I-1) > >2.*3.»GR d > /R 
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GO TO 58 

54 KAPPA2=(GR(2)-0.)/2.+3.*GR(l)/R 
GO TO 58 

56 KAPPA2=GR ( N04M2 ) -GR (N04-3 ) *3.*GR (N04M2)/R 
58 KAPP A2--KAPPA2 

IF<KAPPA2.LT.0.) KAPPA2=0. 

KAPPA ( I ) =SGRT (KAPPA2) 

PRINT 6l*I*R»GR(I) » KAPPA { I > 

60 CONTINUE 

61 FORMATUH *I3*F5.1 .2E16.5) 

PRINT 64 

64 FORMAT ( 1 HO* ANGULARLY AVERAGED QUANTITIES FOR THE END OF THE L'ASt C 
1YCLE WHICH HAVE BEEN COMPUTED AT CELL CENTERS OF MASS AND FOR AyER 

2AGED R.R*) 

PRINT 68 „ - 

68 FORMAT ( 1H * I NO. PART. RNEW MASS DENS 

I SPECIFIC IE PRES ROOT*) 

DO 80 1=1 .N04M2 
RNEWR < I ) =0 • 

MASSR ( I ) =0 • 

PRESR (I)=0. 

GRR=0, 

RDOTR=0. 

NNTH=NTH( I) 

DO 73 J=1 » NNTH 

IF(MASS(I,J) .EQ.O.) GO TO 73 

XNEW=XNWAV ( I*J) 

YNEW=YNWAV(I.J) 

RNEW=SQRT (XNEW*XNEW+YNEW*YNEW) 

RNE'WR ( I ) =RNEWR ( I > .MASS ( I , J> *RNEW 
MASSR ( i ) =MASSR { I ) .MASS ( I , J) 

PRESRU)=PRE5R(I)*mASS( I.J)*IE(I.J> 

X=CENTER+XNEW 

Y=CENTER+YNEW 

IX=X 
I Y=Y • 

XX=IX-X 
YY=I Y-Y 


RAD3=RNEW*RNEW#RNEW 
IF(RNEW.LT.RCC) RAD3=RCC#RCC*RCC 

r*- < YY+ 1 ) * ( ( XX+ 1 • ) {RHO (IX*ltIY)-RH0(IX-l.lY))-XX»CRH0(IX+2.lY> 

iho’uxIiy! > >-yy»uxx*i.>*<rho'1x*i.iy*i>-rhoux-i.iy*i> >-xx*<rho< 
2*2*IY«-1)-RH0<IX»IY*I>) >-GMC*XNEW/RAD3 

r,v- i !tv + l )* { (YY + 1 . )* (RHO( Ix» I Y + l )-RHO ( IX» IY-1 > ) -YY* (RHO ( IX* I Y*2) 
1 rioT IX * I Y ) ) > -XX* ( ( YY+ 1 • ) * (RHO (IX*l*IY*l(-RHO(IX+l«IY-l>)”YYM RHO ( 
2 + l*IY + 2)“RH0(IX+l*IY) ) )-GMC*YNE*/RAQ3 
GRR=GRR+MA5S ( I « J) * (XNEW*GX*YNEW»GY) /RNEW 
RDOTR=RDOTR+MASS(I »J)*RDOT (I.J> 

73 CONTINUE 

IF (MASSR ( I ) • EQ . 0 • ) GO TO 75 
NPARTR=MASSR ( I ) /DMG 


RNE^R ( I ) =RNEWR ( I ) /MASSR { I ) 

DENSR ( I > =MASSR < I ) / (NNTH«AREA ( I > *DT2> 

IER (I)=PREbR{I)/ (MASS r (I)*DT2) 

PRESR { I ) =GM1 *PRESR ( I ) / (NNTH*AREA ( I ) tt DTA) 

GR (I ) =GRR/ (MASSR ( I) *DT2 ) 

RDOTR=RDOTR/ IMASSr ( I )*DT) 
MASSR(I)=MASSR(I)/0T2 


-R 

IX 


-R 

IX 


GO TO 78 
75 NPARTRaO 


I RNEWR(I)-0a 

! densru>=o. 

| IER ( I ) -0 «. 

PRESR ( I ) =0 • 
f GR ( I ) = 0 « 

78 PRINT 81 , 1 ,NPARTR.RNEWR ( I ) , MASSR ( I ) •DENSR ( I ) * IER ( I ) .PRESR < I > * 

I I rdotr 

I 80 CONTINUE 

I 81 FORMAT (1H < 13, I6.6E16.5) 


1 PRINT 83 

\ 83 FORMAT <1 HO* I 

1.THDTO2-THDT02 
i DO 100 I=1*N04M2 

| R2MDR=0. 

R2NWR=0. 

PTHR=0. 

h> 


R2MID 

DP/M 


R2NEW TDOT RNEW 

GRtCNTR MASS) Q (GR (CELL CNTR>>*) 


“» tow vet S syj 


NNTH-NTH(I) 

DO 85 J=1,NNTH 

XF(HASS(I f J).EQ.O.) GO TO 85 
R2MID=R2MDAV(I* J) 

R2NEW=R2NW Av ( I ♦ J ) 

R2MDR=R2MDR +MASS (I*J)*R2MlD 
R2NWR=R2NWR + MA55 ( I » J) *R2NEW 
PTHR=PTHR*HASS ( I * J) »R2MID*TD0T ( I * J) 

85 CONTINUE 

IF (MASSR ( I ) ,EQ. 0 • ) GO TO 90 
R2MDR=R2MDR/ (MASSR ( I ) *DT2 ) 

R2NWR=R2NWR/ (MASSR ( I } *DT2 ) 

TDOTR=PTHR/ (MASSR ( I ) *R2MDR*DT3) 

THDT02=PTHR/ (MASSR ( I ) *R2NWR*DT3) 

RTD0T2=RNEWR ( I ) »THDT02* THDT02 
IFU.EQ.1J GO TO 86 

IF ( I .EQ.N04M2) GO TO 87 
PIP=(PRESR(I) + preSR{I+1) )/2. 

IF(MASSR<I*1).EGU0.5 PIPssO. 

PIM= (PRESr ( i ) ❖PRESR ( I-l ) )/2. 

IF(MASSR(I-I).EO.O.) P iM=0 . 

GO TO 88 

86 PIP= (PRESR ( I ) +PRESR (1+1) )/2. 

IF (MASSR ( I ♦ 1 ) ,eq»o«> pip=o. 

PIM=PRESR(1) 

GO TO 88 

87 PIP=0. 

PlM=(PRE5(I)*PRES(I-i) }/2o 
IF (MASSR (I-1).EQ*0») PIM'-O. 

88 DPOM=(PIM-PIP)/MASSR(I) 

IF(IER( I) .LT.O.) IER ( I ) =0 • 

QR=KAPPA(I)»SQRT(IER(I) )/(3.36*DENSR(I) ) 

GO TO 95 

90 R2MDR=0. 

R2NWR=0. 

TDOTR=0. 

RTDOT2-0. 

DP0M=0. 

QR = 0 ■ 

95 PRINT 101*1 *R2MDR*R2NWR* fDOTR » RTD0T2.DP0M* GR ( I ) *QR 
100 CONTINUE 


101 FORMAtUH »I3»7E16.5) 

120 DO 12? I=1*N04M1 

NNTH=NTH(I) 

DO 122 J=1*NMTH 

IF (MA5S(I*J) • EQ . 0 • ) GO TO 121 

PRES { I « J ) =GHl tt I£ (I*J)*MAS5{I*J)/AREA(I) 

IFUEi I* J> .LT.O.) PRES ( I , J) -0 • 

GO TO 122 

121 PRES<I,J>=0. 

122 CONTINUE 
iete^t=o 
nhivel=o 
5UMPEF=0 • 

sumpes=o. 

DO 200 I“1 *N0AM2 
ARCIP=I*DTH(I) 

ARCIMs(I-l)*DTH(I> 

IPTESt=0 

IF(I.EQ.l) IPTEST=1 
IF ( I «EQ*2) IPTEST=1 
IFtI.EQ.4) IPTEST=1 
IF(I.EQ.B) IPTEST=1 
IF(I.EQ.16) IPTEST=1 
IMtESt=0 

IFU.EQ.2) IHTEST=1 
IFU.EQ.3) IMTEST^l 
IF(I.EQ.5> IMTEST^l 
IF ( I • EQ *9 ) IMTEST = 1 
IF ( I«EQ* 17) IHTESt=1 
NNTH=NTH{I) 

DO 200 J-l *NNTH 

IF (MASS ( I * J) , EQ • 0 . ) GO TO 200 

XNEW=XNWAV ( I * J) 

YNEW=YNWAV<I»J) 

RNEW=SQRT iXNEW^XNEW+YNEW*YNEW> 

x=center+xnew 

y=center+ynew 

IX=X 



IY = Y 

XX=IX-X 

YY=IY-Y 

RAD3=RNEW*RNEW*RNEW 
IF(RNEW.LT.RCC) RAD3=RCORCC*RCC 

SX=(YY+1.)*( (XX+1.)*(RH0( IX*l*IY)-RHO<lX-ltIY> )-XX*(RHO< IX*2* IY)-R 
1H0(IX,IY) ) )-YY*( <XX+l.>*(RHO(IX*l*IY>l)-RHO(IX-lflY+l) >-XXMRHO{IX 
2 + 2*IY+l)-RH0<IXtIY + l) > ) -GMC*X\IEW/RAD3 
GY=(XX*1.)*( <YY+l.>*(RHO(lX*IY+l)-RHO(lXtIY-l> > -YY* <RHO ( IX t I Y*2 ) -R 
1H0(IX.IY)))-XX*( (YY*1.)*(RH0(IX+1*IY*1)-RH0(IX+1*IY-1) >-YY*<RHQ{IX 
2+l*IY+2)~RH0(IX+lfIY) ) ) -GMC*YMEW/RAD3 
DX=-XX 
DY=-YY 

D11=1-DY-DX*DX»DY 
D12=DY* (1-OX) 

D21=Dx a ( 1-DY) 

D22”Dx*DY 

SUMPESaSUMPES*MASS C I » J) * (Dll^RHO (IX, IY) ♦D12*RH0 ( IX* IY + 1 ) 

1 ♦D2l#RH0<lX*l*IY>*D22 tt RH0(IX+l*IY+l> ) 

EGC=GMCP1-GMCP2*RNEW*RNEW 
IFtRNEW.GT.RCC) EGC=GMC/RNF.W 
SUMPEF=SUMPEF*MASS { I , J ) »EGC 
IF(IPTESt.EG.I) GO to 130 
PIPa(PRES(I, J>*PRES(I*lfJ) )/2. 

IF(MASS(I+1, J) .EQ.O.) PIP=0. 

DEIP=-ARCIP*PIPMRDOT< I*J) ♦R00T<I*1*J> )/2. 

GO TO 1 AO 

130 PIPl“{PRES(If J) +PRES(I+1,2*J-1) )/2. 

IF (MASS { 1+1 ♦ 2*J- 1 ) .EQ. 0 • ) PIP 1=0. 

PIP2=(PRES( I , J) *PRES{ 1+1 ,2»J) )/2. 

IF (MASS { 1 + 1 t2*J) .EQ.O. ) PIP2=0. 

PIP=tPIPl+PIP2)/2. 

0EIP=-.5*ARCIP*PlPl tt tROOjdt J) *RDOT( I + 1*2*J-1) )/2. 

1 -.5*ARCIP # PIP2 tt (RDOT<IfJ) ♦RDOT (Id t2*J) )/2. 

140 XFCI.NE.l) GO TO 145 

JP2=J+2 

IF { JP2.GT .4) JP2=JP2-4 
PIM=(PRES{lfJ> +PRES<1,JP2) )/2. 

IF(MA5S(lfjP2) .EQ.O.) PIH=0. 

DEIH=0. 


145 


160 


JJ=J 

IF t IMTEST.EO.l ) JJ- ( J* 1 ) /2 
PIM= (PRES( I ,J)*PRE6(I-l*JJ))/2. 

IF(MA5S(1-1, JJ) .EQ.O.) PIM=0. 

DEIM = ARCIM fr PlM* (ROOT ( I * J) +RDOT ( 1-1 * JJ> > /2. 
JP1=J*1 

IFIJP1.GT.NNTH) JP1*1 

PJP= (PRES ( I , J) +PRES ( I * JP1 ) ) /2. 

IF (MASS( I,JP1) .EQ.O.) PJP = 0. 

QEjp=-,pjp» ( I -.5) * (TDOT IItJ)*TDOT(I*JPl))/2. 
JH 1 = J- 1 

IF(JMl.LT.l) JM1-NNTH 
PJM*(PRES(I,J)*PRES(I.JM1) )/2. 

IF(MA5S(I*JM1).EQ.0.) pJm = °* ,,, mllw o 

DEJM=PJM* ( I-.5) * (TDOT < I » J) ♦TOOT (I*JM1) >/2. 

GR= (XNEW*GX+VM£W*GY) /KNEW 

FR=GR+ (P IM-P IP ) /MASS ( I * J > 

GTH= (XNEW*GY-YNEW*GX) /RNEW 

FTH=GTH+ (PJM-PJP)/«ASS< I* J) 
R2MID=R2MDAV(I#J> 

R2NEW=R2NWAV(I*J) 

PTH=R2MID»TD0T(I*J) 

TORQ=RNEW*FTH 

pthot=pth+torq 

THDT02=.5*<PTH«-PTHDT> /R2NEW 
RDDOT 1 1 , J) =FR*RNEW*THDT02 W THDT02 
RDOT0t=RDOT n ♦ J) +RDDOT ( I ♦ J) 

RD0T34=RD0T ( I * J) ♦ .75*RDD0T ( 1 1 J> 

R20T=R2NEW+RNEW*R00T34+ .25*RDDT34*RD0T34 

IF (R2DT .EQ.O » > R2DT=1. 

TD0TD T =PTHDT/R2DT 


TF-fT n — ifci,jj + .5* (ROOT ( I , J) *RDOT ( I * J) ♦R2MID*TDUT t 

1 IE a -RD0TDT*RD0TDT-R2DT*TD0T0T*TD0TDT)+GR*(RD0Ta*J 
2 +GTH*RNE^*THDT02+ { DElP + DEIM + DE JP + DEJM > /MASS ( I * J 

IF ( IE ( I * J) *LT .0 • ) IETESjssl 
IFtMASS(ItJ) »LT.5.0*DMG) GO TO 1B5 
ABROOT=ABS(ROOTDT) 


J) ♦R2MID*TD0T ( I « J) *TD0T < I » J* 
> +ROOTDT) /2* 


ABTDOT=ABS{TDOTDT) 

IF <ABRDOT.GE..375.0R.ABTDOT.GE..375*DTH<I ) ) NHl VEL=NHI VEL + 1 

IF(ABRDOT c LT.l .0. AND. ABTOOT.LT. DTH(I) ) GO TO 185 

CY=CYY 

IPL0TG=1 

IPRNTG=1 

PRINT 165»RD0TDT*T00TDTv v * J 

165 FORMAT(1HO*RDOT=*E16.8* AND TD0T=*E16.8* FOR I-*I3* AND J=»I3) 
165 PTHET A ( I * J> =TDOTDT 
200 CONTINUE 

IF{NHIVEL.GT.O> PRINT 205.NHIVEL 

205 FORMAT { 1H0*NUMB£R OF CELLS WITH ,GT. 5 PART I CLE5 AND WITH SCALED V 
1EL0CITY COMPONENTS .GT. .375. (CELL DIMENSION) = *15) 

IF(IEtESt.EQ.1.AND.CYY-CyY/ 5*5.E3-.0> CALL’ NEGIE 
PTH=SUMPTH/DT3 
TM=SUMMAS/DT2 
INE=SUMIE/DTA 
KEsSUMKE/DT4 
PESNEW=~SUMPES/DT4 
PEFNEW=-SUMPEF/DT4 
IF(CYY.EQ.l) PESOLD=PESNEW 
IF(CYY.EQ.l) PEFOLD^PEFNEW 
PES= (PESOLD+PESNEW) /2. 

PEF =(PEF0LD+PEFNEW>/2. 

pesold=pesnew 

PEF0LD=PEFNEW 

PE=PES-*-PEF 

te=pe*ine+ke 

PRINT 232 

232 FORMAT (1H0 *THE VALUES ON THE NEXT TWO LINES ARE FOR THE END OF TH 
IE LASj CYCLE*) 

PRINT 234»PES,PEF,PTH,TM 

234 FORMAT < 1 H *PES=*E 14 . 7* PEF=*E14.7* PTH=»E14.7* TM=*E14.7> 

PRINT 236, PE, INE,K£,TE 

236 FORMATUH *PE=*E14.7* IE=*E14.7* KE=*E14.7» TE=*E14.7) 

ICYPLT=CYY-CYMIN 
SAVPTH(ICYPLT)=PTH 
SAVPE ( ICYPLT ) =PE 


rt 

i 




SavIE(ICYPLT)=INE 
SAVKE ( ICYPLT ) =KE 
SAVTE(ICYPLT>=TE 
DO 245 I=1*N04M1 
MNTH=NTH(I ) 

DO 245 J=1 »NMTH 

IF (MASS(IfJ) . EQ« 0 , ) GO TO 240 

TDOT (I»J)=PTHETA(I*J> 

240 PRAD ( I » J) =0 » 

PTHET A ( I , J ) =0 • 

etotal ( I ♦ J) =0 • 

MASS ( I , J ) =0 • 

XNwAV ( I * J ) -0 • 

YNWAV(I*J)=0. 

R2MDAV (I * J ) =0 * 

R2NWAV ( I ♦ J) “0 * 

245 CONTINUE 

DO 248 1=1 *N02 
DO 248 J=l*N02 
RH0( I , J)=0. 

248 CONTINUE 
PRINT 250 

250 FORMAT <1H0*LAST PARTICLE OF 1ST 10 SETS OF NBSG GAS PARTICLES*) 
PRINT 260 

260 FORMAT (1H ^NUMBER X-CENtER Y-CENTER R 

1 THETA ROOT TDOT*) 

NOUT=0 
NPART=1 
NS2 = 0 

270 READ<JTG) RADIUS, THETA 
DO 350 IS=1,NBSG 
ROLD=RADIUS(IS) 

IF(ROLD.GT.RMAX) GO TO 340 

IROLD=ROLD*l 

TOLD=THETA(lS) 

JTOLD=TOLO/DTH(IROLD)^1 

RD0T34=RD0T (IROLD, JTOLD) ♦.75*RDD0T ( IROLD, JTOLD) 

RMID=ROLD+ *5*RD0T34 






RDOTT=RDOT (I ROLD, JTOLD ) +RDDOT ( IRQLD, JTOLD) 
rnew=rold*rdott 

IF (RNEW.GT.RMAX) GO TO 340 
TDOTT=TDOT ( IROLO* JTOLD) 

tnew=told*tdott 
IF(RNEW) 273 *274*275 

273 RNEW=-RNEW 
tnew=tnew+pi 
GO TO 275 

274 RNEW=1.0E-8 

275 IRNEW=RNEW*1 

IF ( TNEW *LT • 0 • ) TNEW=TNEW + TWOPI 

if ctnew.ge.twopi ) t new=tnew-twopi 

JTNEW=TNEW/DTH(IRNEW) *1 

DPR=DMG*RDOTT 

R2MID-RMID*RMID 


DPTH=DMG*R2MID*TDOTT 

DE=DMG <J ( IE ( I ROLD* JTOLD) + .5MRD0TT*RD0TT + R2MID*TD0TT*TD0TT 

PR ADC IRNEW, JTNEW )=PR ADC IRNEW, JTNEW) +DPR 

PTHETA ( IRNEW, JTNEW) =PTHETA ( IRNEW * JTNEW ) +OPTH 

ETOTAL ( IRNEW, JTNEW )=ET0TAL< IRNEW, JTNEW) ♦DE 

MASS ( IRNEW, JTNEW) =MASS ( IRNEW, JTNEW) +DMG 

R2MDAV < IRNEW, JTNEW) =R2MDAV ( IRNEW, JTNEW) *R2M ID 

R2NWAV ( IRNEW, JTNEW )=R2NWAV ( IRNEW, JTNEW ) ♦RNEW^RNEW 

XNEW=RNEW*COS(TNEW) 

YNEW=RNEW*SIN (TNEW) 


IXNEW=CENTER+XNEW+.5 

jynew=center + ynew+.s 

RHOCIXNEW, JYNEW)=RHO(IXNEW, JYNEW ) ♦ ,.5*DMG 


XNwAVtlRNEW, JTNEW ) 
YNWAVCIRNEW, JTNEW) 


=XNWAV ( IRNEW, JTNEW) +XNEW 
=YNWAV < IRNEW, JTNEW ) +YNEW 


GO TO 345 

340 NOUT=NOUT+1 
RNEW=999. 


) 


TNEW=999. 

XNEW=999. 

YNEW-999. 

TD0TT=999* 

RD0TT=999. 




3A5 RADIUS(IS)=RNEW 


THETA ( IS) =TNEW 

350 CONTINUE 

WRITE(JSG) radius, theta 

IF (NPART.LE. 10)PRINT 360 , NPART » XNEW , YNEW »RNEW , TNEW , RDQTT ♦ TDOTT 


360 FORMAT < 1H ,I6,6E!6.8) 

NPART =NPART + 1 
IFIIPLOTG.EO.O} GO TO 370 
DO 365 IS=1,NBSG 
R=RADIUS(IS) 

TH=THETA(IS) 

XPLOT (IS)=R*C0S(TH) ♦CENTER 
YPLOT (IS)~R*SIN(TH) +CENTER 


365 CONTINUE 


Q=0 

IF (NS2.EQ.NMKG) Q=1 

CALL DDIPLT (Q*ING«NBSG*XPLOT»YPLOT»XMING*xMAXGtYMING*YMAXG» 
I 2,XPG»1*YPG»13*ITAPXG> 

370 NS2=NS2+NBSG 

IFtNS2.LT.NBRG) GO TO 270 

REWIND JTG 
REWIND JSG 
JSAVE=JSG 
JSG=JtG 

jtg=jsave 
PRINT 373*N0UT 

373 FORMAT {1H0*NUMBER OF PARTICLES 0UT5IDE OF MESH=*I6) 

summas=o. 

sumpth=o. 

sumke=o. 

SUMIE=0 • 

DO 3B0 I=I*N04M1 
NNTH=NTH( I ) 

DO 380 J=1»NNTH 

IF (MASS ( I ♦ J) . EQ.O . ) GO TO 376 

NUM8ER=MASS { I , J) /DMG 

XNWAV ( I , J> =XNWAV ( I, J) /NUMBER 

YNWAV (I»J)-YNWAVd* J> /NUMBER 

R2MDAV ( I * J) =R2MDAV ( I » J) /NUMBER 


R2NWAV ( I * J) =R2NWAV < I * J) /NUMBER 
SUMMAS=SUMMA5>MASS (I, J) 

; RDOT(I.J)=PRAD(I*J)/MASS<I,J) 

I SUMP TH=SUMPTH+PT MET A { I ♦ J ) 

= IF(R2MDAV(I,J>.EQ.0.) R2MDAV(I*J>=1.0E-8 

; TD0T(I.J)=PTHETA(I»J)/(MASS<I,JJttR2MDAV(I*J) ) 

I «E=.5*MASS(I,J)*(RDOT(I*J)*RDOT<I»J> 

1 +R2MDAv(I*J>*TDOT<I*J>#TDOT<I*J) ) 

I 5UMKE=SUMKE+KE 

IE(I*J)=ETOTAL(I»J)-KE 
( SUmIE=SUMI£*IE ( I* J) 

| IE(I*J)=IE(I*J)/MASS(I,J) 

I; GO TO 380 

376 RDOT ( I ♦ J } =0 * 

\ TDOT < I » J ) =0 * 

I IE ( I * J)=0» 

380 COMTINUE 

WR ITE (8 ) RDOTtTDOT. lE,MASS»XNrtAV*YNWAV*R?MDAV*R2NWAV 
j REWIND 8 

IF ( IPRNTG. NE* 1 ) GO TO 485 
PRINT 390 

I 390 FORMAT (1H0»J» (RDOT ( I ♦ J) , 1=1 . 15) *) 

400 FORMAT ( 1H ,I4»15F8.4> 

I DO 410 J=l,64 

PRINT 400* J* (RDOT ( I * J) * 1 = 1 ♦ 15) 

I: 410 CONTINUE 

I PRINT 420 

420 FORMAT (1H0»J. (TOOT ( I * J) ♦ 1 = 1 ♦ 15) 

1 00 430 J=l*64 

PRINT 400, J, (TD0T(ItJ)»I=l»15) 

% 430 CONTINUE 

| PRINT 440 

l 440 FORMAT (1H0*J, (IE ( I , J) , 1=1 , 15) ») 

| 450 FORMAT (1H ,I4,15E8.1) 

I DO 460 J=1 , 64 

PRINT 450»J, (IE(I,J) '*1 = 1*15) 

\ 460 CONTINUE 

f. PRINT 470 

| 470 FORMAT (1H0»J, (MASS( I , J) *I=:l*15) tt > 

% 

S> 

% 


DO 480 0=1*64 

PRINT 450 *J* ( MASS ( I * J ) *1-1*15) 


480 CONTINUE 

485 IFUPLOTG.NE.l ) GO TO 600 
C PLOT VELOCITY VECTORS 
XN02=N02 
A=0. 

CALL DDIPLT (0* ING*1 *A* A*0*XN02*0*XN02*2*XPG*1*YPG» 13*ITAPXG) 

XYSCAL=10./N02 
VSCALE=.001/DT 
DO 500 I=1*N04M2 
NNTH=NTH(I) 

DO 500 J=I*NNTH 

IFtMASS (I, J) .EG.O.) GO TO 500 


XNEW=XNWAV(I*J) 

YNEW=YNWAv<I*J) 

RNEW=SQKT (XNEW»XNEW+YNEW«YNEW) 
VX=-YNEW*TDOT(I*J>*XNErt*RDOTU*J)/RNEW 

VY=XNEW*TDOT ( I * J) +YNErf*RDOT ( I * J) /RNErt 

XM lD=XNEW- , 5 tt VX 

YMI0=YNEW-.5*VY 

X=XMID*CENTER 

y=ymid+center 

XA=XYSCAL*X 

YA=XYSCAL*Y 

XB=XYSCAL*(X+VSCALE*VX) 

Y8=XY5CAL* (Y+VSCALE*VY) 
call PARR0W(XA*YA,XB*YB*1 ) 

500 CONTINUE 

CALL DDIPLT u * ING.lf A*A»0*XN02*0*XN02*2*XPG*1*YPG*I3 * ItAPxG) 
C HAKE A CONTOUR PLOT OF THE DENSITY 

CALL DDIPLT (0* ING*5*XPR*YPR*XHIN4*XMAX4*YMIN4*YHAX4*3*XYL'{ I ) 
1 2,PLH(1> *14*ITAPXG) 

DX = 0. 

DY-1 . 

DO 515 J=2.N02 
K=0 

DX=DX*RCOS 

DY=DY+RSIN 

DU=0. 

TV DV=1. 

J . 


***'*'■■ 



00 510 l=2tN02 
K=K ♦ 1 

DU = DU4-RC0S 
DV=OV*RSIN 
XC=I-CENTER 

YC=J-CENTLR 

IRsSQRT (XC*XC*YC*YC> ♦ ! 

IF ( IR.GT .NQ4M2) GO TO 506 

TH=ATAN2(XC,YC) 

IFdH.LT.O.) TH=TH*TW0PI 

JT1 = TH/DTHUR) +1 
ZI=MAS5(IR,JT1)/AREA(IR) 

TH=ATAN2 (YCfXC) 

IFtTH.LT.O.) TH=TH*TW0PI 
JT2=TH/DTH(IR) + 1 
Z2=MASS(IR,JT2)/AREA(IR> 

GO TO 508 
506 Z1 = 0 . 

Z2=0. 

508 XPL0T1 (K ) = J«-DU 

YPL0T1 (K)=SPHIM»Z1*DV 
XPL0T2 (K) =I+DX 
510 YPL0T2 (K) =SPHIH*Z2+DY 

CALL DDIPLT(0,ING,K,XPLOt2,YPLOT2»XMINA*XMAXA»YMINA*YMAX4, 

1 3,XYL(1) »2*PLM{1) *14,lTAPXG) 

CALL DDIPLT (Of ING.K»XPLOH »YPL0T1 »XMlN4»XMAX4*YMIN4*YMAX4t 
1 3,XYL(1) f2fPLH(l) f lAflTAPXG) 

515 CONTINUE 
K-0 

DO 520 I=2fN02 
XI = I 
K=K + 1 

XPL0T2(K)=XI*RC0S 

YPL0T2(K)=1.*RSIN 

520 CALL^DDIPLT (If iNGfKf XPLOTSf YPLOTHf XMINAf XMAXAt YMINAfYHAXAf 

1 3,XYL(1> *2fPLM(l> t 13,ITAPXG) 

C MAKE A CONTOUR PLOT OF THE TEMPERATURE 

CALL DDIPLT (Of INGt 5f XPRfYPRf XMINAf XMAXAfYMINAf YMAX4t3f XYL'(l) « 

1 2,PLT(1) f 14f ITAPXG) 

DX=Q. 



DY-1. 

DO 540 J=2 *N02 
K-0 

DX=DX*RCOS 
DY=DY*f:3lN 
DU = 0 . 

DV = 1 • 

DO 530 1=2 *N02 


X=K + 1 

DU=DU*RCOS 

DV=DV+R51N 

XC=I-CENTER 

YC=J-CENTER 

IR=SQRj (XC*XC*YC*YC) +1 

IF ( IR.GT.N04M2) GO TO 524 

TH=ATAN2(XC»YC) 

IFITH.LT.O.) TH=TH*TW0PI 

JT1=TH/DTH(IR)+1 

Z1=IE ( IR* JT1 ) 

TH=ATAN2 (YC»XC> 
IFCTH.LT.O.) TH=TH+TWOPI 
JT2 = TH/DTHUR) *1 
Z2=IE{IR*JT2) 

GO TO 526 
524 Z1=0. 


>3 

| 

.'if' 

::i: ; 

I 

ft 


Z2 = 0 « 

526 XPLOT1 (K) =J+DU 

YPLOT1 (K)=SpHIT*Zl*DV 


XPL0T2(K)=I*DX 


YPLOT2 (K) =5PHIT*Z2+DY 

CALL DDIPLT(0,ING,K,XPLOT2,YPLOT2,XMIN4, 


uti i u / UltTlI/. 


1 3»XYL<1) »2«PLT(1) »14*ITAPXG) 

CALL DOIPLT (Of INGtK*XPLOTl *YPL0T1*XMIN4*XMAX4-*YMIN4*YMAX4« 

1 3,XYL(1) »2.PLTU> »14,ITAPXG) 


540 CONTINUE 


l <=0 ■ 

DO 545 I=2*N02 
XI = I 


K=K>1 

XPLOT2(K)=Xl+RCOS 








YPL0T2<K>=1.+RSIN 
545 CONTINUE 

CALL DDIPLTU ♦ING,K,XPLOt2,YPLOT2»XMIN4»XMAX4*YMIN4*YMAX4« 

1 3»XYL(1) f2*PLT(l) »13 ,ItAPXG) 

C MAKE A CONTOUR PLOT OF THE PRESSURE 

CALL DDIPLT <0» ING»5fXPR»YPR»XMIN4,XMAX4*YMIN4,YMAX49 3»XYL (l) ♦ 
1 2 *PLP (1 ) *14*ITAPXG) 

DX=0. 

DY = 1. 

DO 560 J=2 *N02 
K = 0 

DX=DX*RCOS 

DY=DY+RSIN 

DU=0. 

DV- 1 * 

DO 550 1-2 *N02 
K-K* 1 

DU=DU>RCOS 
DV-DV+RSIN 
XC=I-CENTER 
YC=J-CENTER 
IR=SQRt<xC»XC*YC*YC> ♦! 

IF ( IR.GT.N04M2) GO TO 54? 

TH=ATAN2(XC < YC) 

IF (TH.LT .0 » ) TH=TH*TWOPI 
JT1=TH/DTH( IR) +1 

Z1=IE<IR,JT1)*MASS<IR,JT1)/AREA<IR) 

TH=ATAN2(YC.XC) 

IF(TH.LT.O.) th=th*twopi 
JT2=TH/DTH(IR)*1 

Z2=IE(IR,JT2)*MASS<IR,JT2)/AREA(IR) 

GO TO 549 
547 Z1 = 0. 

Z2=0. 

549 XPLOT1 (K)-J+DU 

YPLOT1 (K)=SPHIP»Z1*DV 


XPLQT2 ( K ) -I +DX 
550 YPLOT2(K>=SPHIP*Z2*DY 

CALL DD IPLT (0* ING » K , XPLOf 2»YPL0 t2, XM IN4, XM AX4 » YM IN4 » YM AX4 » 
1 3,XYL(l)t2»PLP(l)flA,lTAPXG) 

CALL DDIPLUO* ING*Kt XPLOTl , YPLOT 1 , XM I N4, XM AX4, YM IN4 , YMAX4 » 
1 3,XYL(l)*2»PLP<l>*14,lTAPXG) 

560 CONTINUE 


K=0 

DO 570 I~2 »N02 
X 1 = I 


K=K ♦ 1 

XPLOT2 (K)=XI+RCOS 
YPLOT2 (K)=I.*RSIN 


570 CONTINUE 

CALL DDIPLT(l*INGfKtXPL0'T2,YPL0T2*XMlN4fXMAX<nYMIN4*YMAX4* 

1 3 f XYL(l ) »2»PLP U ) *13t ITAPXG) 

600 IF(CYY.LT.CY) RETURN 

C IF(CYY.GE.CY) WRITE INFORMATION TO BE SAVED ON TAPEI 
NS2-0 

610 REAO(JTG) RADIUS, THETA 
WRITE ( 1 ) RADIUS, THETA 
NS2=N52«-NBSG 


IFtNS2.LT.NBRG) GO TO 610 

WRITE (1 ) RDOT » TOOT » IE, MASS, XNWAV»YNWAV»R2MDAV»R2NWAV»SUMMAS, 

1 SuMPTH,SUMIE,SUMKE,PES0LD,PEFOLD,CYY,DT*DT2,DT3,DT4,DMG,GMC, 

2 GMCP1 .GMCP2»SPH I M,SPh I T»S PH IP, PLM (2> *PLT (2) »PLP (2) ,EMAX,EMIN» 


2 GMCP1 

3 PTHMAX,PTHMIN,SPHI,PL (2) 
REWIND 1 


RETURN 


end 


-jo ui u ro 


SUBROUTINE negie 

COMMON/ AILC0M/I2A, I TESt»N«N02»CY, CYY»RH0 (32,32) 

COMMON/GA5COM/NTH( 15) «DTH(I5) ,AREA(15) ,NPLOTG,NPRNTG,NRING,N04» 

1 N04M1 »N04M2»RCC,GMC»GMCP1»GMCP2,DMG*Dt,DT2»DT3»DT4,NRTCEL» 

rmax, center, SUMPTH,SUMM AS, SUMIE,SUMKE,PES0LD,PEF0L0*XMG, 

SAVPTH (200) , SAVPE (200) ,SAVIE (200) ,SAVKE{200) ,SAVTE(200) , 
MASK1,N8SG,NBRG, JTG, JSG, IMG(2) ,XMING,XMAXG,YMING,YMAXG,X»G(2) , 
YPG,ITAPXG,GM1,NMKG,XYL(3) ,X?R(5) ,YPR(5) *XMIN4,XMAX4, YMIN4, 
YMAX4,PLM(2) ,RC0S,RSin,SPhIM,PLT(2) , SPHI T ,PLP ( 2 ) »SPHIP»TW0PI » 
PTHMIN,PTHMAX,EMIN»EMAX,CYMIN,PI,PL (2) ,SPHI,MTESt,ITAPE 
COMMON/NEGCOM/IE(15,64) ,MASS(15,64) ,IEC0R(15,64) 

REAL IE,IEC0R,IE1, IE2,IE3,IE4, IE5, IE6, IE7, IE8, IE9,MASS 
DO 5 I=1,N04M1 

NNTH=NTH(I) 

DO 5 J=1,NNTH 
IECOR«I,J>=0. 

IE(I,J)=MASS(I,J)*IE(I,J) 

5 CONTINUE 
IE ( 1 , 8 ) =0 . 

IE<1»9)=0. 

IECOR ( 1 , 8 ) =0 • 

IECOR ( I , 9) =0 • 

DO 80 1=1 ,N04M2 
NNTH=NTH(I) 

IRTE5 T =1 

IF(I.EQ.l) IRtESt= 2 
IFd.EQ.2) IRTESt= 3 
IFd.EQ.4) IRTEST=4 
IF(I.EQ.8> IRTESt = 4 
IF (I «EQ p 16) IRtESt=4 
IF(I.EQ.3) IRTEST=5 
IF (IpEQp5) IRTESt=5 

IFd.EQ.9) IRTESt=5 
IF ( I *EQ«17) IRTESt=5 
DO 80 J=1,NNTH 



IF ( IE ( I * J) .GE.O . ) GO TO 80 
GO TO (10*20*30,40.50) IRTESY 

10 11 = 1+1 
J1=J 

12 = 1+1 
J2=J+1 
13=1 
J3=J+ 1 
14=1-1 
J4= J+ 1 
15=1-1 
J5=J 
16=1-1 
J6=J-i 
17=1 
J7= J- I 
18=1+1 

J8= J- 1 

19=1 

J9=9 

IF(J.NE.l) GO TO 15 
J6=NNTH 

J7=NNTH 
J8=NNTH 
GO TO 60 

15 IF ( J.NE.NNTH) GO TO 60 
J2=l 
J3= 1 
J4=l 

GO TO 60 
20 11=1 
J1=J+1 
12=1 
J2=J+2 
13=1 

J3=J+3 

14=2 


J4=2*J-2 

15=2 

J5=2* J- 1 
16=2 
J6=2*J 
17=2 

J7=2* J* 1 
18 = 1 
J8=8 
19=1 
J9=9 

IF ( J1 ,GT*4) Jl=Jl-4 
IF(J2.GT.4> J2=J2-4 
IF<J3.GT.4) J3=J3-4 
IF(J7,GT.8) J7=J7-8 
IF<J4.LT.l) J4=J4+8 
GO TO 60 
30 11=3 
J1=2»J 
12=3 

J2=2* J«- 1 
13=2 

J3=J+ 1 
14=1 

J4=J/2+ 1 
15 = 1 

J5=J/2 

16=2 

J6=J-1 

17=3 

J7=2*J-2 

18=3 

J8=2*J-1 

19=1 

J9=9 

IF(J.NE.l) GO TO 35 
J5=4 




J6=8 
J7= 1 6 
GO TO 60 

35 IFU.NE.8) GO TO 60 
J2= 1 
J3=l 
J4=l 

GO TO 60 

40 11=1+1 
J1=2*J 
12 = 1+1 
J2=2*J4 1 
13=1 
J3=J+ 1 

14=1-1 ' 

J4= J+ 1 
15=1-1 
J5 = J 
16=1-1 
J6=J- 1 
17=1 
J7=J- l 
18=1+1 
J8=2* J-2 
19=1+1 
J9=2*J-1 

IFU.NE.l) 60 TO 45 
J6=NTH(I6) 

J7=NTH (17) : 

J8=NTH(I8) 

. GO TO 60 

45 IF { J.NE.NNTH) GO TO 60 

J2=l 

J3=l 

J4= 1 

GO TO 60 






U'-d 


50 11=1+1 
J1=J 
12 = 1+1 
J2=J+1 
13=1 
J3=J+ 1 
14=1-1 
J4=J/2+ 1 
15=1-1 
J5=J/2 
16=1 
36=3-1 
17=1+1 
J7=J-1 
18 = 1 
J8=8 
19=1 
J9=9 

IF(J.NE.l) GO TO 55 
J5=NTH(I5) 

J6=NTH (16) 

J7=NTH (17) 

GO TO 60 

55 IF< J.NE.NNTH) GO TO 60 
J2=l 
J3= 1 
J4= 1 

60 IE1=IE<I1*J1> 

IF(IE1,LT.0.) IE1=0. 
IE2=IE ( 12* J2) 

IF ( IE2.LT *0 * ) IE2=0. 
IE3=IE ( 13* J3) 

IF (IE3.LT *0#) IE3=0 • 
IE4=IE ( 14* J4) 

IF ( IE4.LT ) IE4=0. 
IE5=IE (15*35) 



IFUE5.LT. 0.) IE5 = 0 • 
IE6 = IEU6,J6> 

IF ( IE6.LT .0. ) IE6=0. 
IE7 = IE { 17. J7) 
IFUE7.LT. 0.) IE7=0 . 
IE8=IE(I8,J8) 

IF ( IE8.LT .0 • ) IE8=0. 


IE9=IE(I9»J9> 


IF < IE9.LT .0. ) IE9-0. 
SUM=IE1 + IE2UE3+IE4* 


IE5+IE6+IE7UE8+IE9 


IF<SUm.LT.-IE(I*J) ) GO TO 75 


QUOTNT = IEU.J)/SUM 

IECORdlf J1) = IEC0R(I1.J1) ♦QU0TNT*IE1 


IECORU2* J2) = IEC0R(I2»J2) ♦0U0TNT*IE2 


IECORU3»J3)=IECOR(I3*J3) +QU0TNT tt IE3 
IEC0RU4, J4) =IECORU^» J4) +QU0TNT*IE4 


IEC0R( I5»J5)=IEC0R(I5, J5) ♦QUOTNT^lES 


IEC0R<I6,J6>=UEC0RU6*J6) ♦QU0TNT*IE6 
IE COR (I7*J7) = I£C0RU7»J7) ♦QUOTNT*IE7 


IEC0R(I8,J8)=IEC0R(I8,J8) ♦QU0TNT*IE8 
IECOR U9*J9) = IEC0R(I9,J9) ♦QUCTNT*IE9 
75 IECOR (I . J) =IECOR (I*J)-IE(I»J) 

80 CONTINUE 

DO 100 I~1 *N04M2 

NNTH=NTH<I) 

DO 100 J=1»NNTH 
IE (MASS (I.J).EO.O.) GO TO 100 
IE(I*J) = UEU»J)*I ECOR ( I , J > ) /MASS (I , J ) 
IF (IE(I*J) .LT-O.) IE < I * J) =0 « 

IECOR < I * J) =0 «, 

100 CONTINUE 
RETURN 
END 




OVERLAY (IFILE«4*0) 

PROGRAM GASPLT 

COMMON/GASCOM/NTH (15) *DTH< 15) , AREA (15) ♦NPLOTG.NPRNTG,NRING*N04, 

1 NC4M1 ♦N04M2 * RCC » GMC ♦ GMCP1 * GMCP2 » DMG* DT tOT2*DT3* DT4* NRTCEL'* 

2 RMAX ♦CENTER*SUMPTH, SUMMA5,SU^ I E»SUMKE,PEGOLDf REFOLD »XMG* 

3 SA VPTH ( 200 ) *5AVPE(200) ♦SAVIE(200) »5AVKE(200) *SAVTE(200) » 

4 MASK1 »NBSG»NBRG» JTG» jSG* IMG (2) ,XMING*XMAXG*YMING»YMAXG*XPG (2) * 

5 YPG. ITAPXG»GM1,NMKG*XYL (3) »XPR<5) » YPR (5) *XM IN4*XMAX4» YMIN4, 

6 YMAX4»PLM(2) »RC0S,rSin,SPHIM*PLT<2) » SphIJ *PLP < 2) .SPHlPtTtfOPI ♦ 

7 PTHMIN*PTHMAX»EMlN*EMAX *CYMIN»PI *PL (2) »SPHl *MTESr * ITAPE 
COM MON/ALL COM/ I 2A* ITESt *N » N02 . CY » CYY « RHO ( 32 » 32 ) 

DIMENSION XDATA(200> *LABLPTH(2) *LABLEN(8) 

DIMENSION ENDCYY (2) *ENDPTH (2) ♦ ENDPE (2) ,ENDIE(2) *ENDKE ( 2 ) * ENDTE ( 2 ) 
INTEGER CY 

CALL PSEUDO * 

CYMAX-CY 

NCYY=CY-CYMIN 
DO 20 Isl,NCYY 
20 XDATA{I)=CYMIN*I 
lablcyy=iohcycles 
LABLPTH(1)=10HANGULAR MO 

lablpth(2)=iohmentum 

LABLEN(1 >=10HENERGY»**C 

LABLEN(2)=10HIRCLE-P0TE 

LABLEN(3)=10HNTIAL»SQUA 

LABLEN(4)wlOHRE~ INTERNA 

LABLEN(5)-10HL»DIAMOND- 

LABLEN(6)=10HKINETIC*TR 

LABLEN (7) =10HI ANGLE-TOT 

LABLEN(8)=10HAL 

ENDCYY U )=CYMIN+1 

ENDCYY (2 ) =CY 

ENDPTH ( 1 ) =SAVPTH ( I ) 

£NDPTH(2)=SAVPTH(NCYY> 



ENDPEC1 )=SAVPE(1> 
ENDPE ( 2 ) =savpe (NCYY ) 
EN0IEU)=SAVIE(1> 

ENDIE(2)=SAVIE(NCYY) 


ENDKE(1 )=SAVKE(1) 


ENDKE(2»=5AVKE(NCYY) 

EN0TE(1)=SavTE(1> 


ENDTE < 2 ) =SAVTE (NCYY > 


IF (NCYY «LT »2) GO TO 30 

CALL DDIPLT <0, ING, 2, ENDCYY,ENDPTH,CYMIN, 


CYMAX , PTHM IN, PThMAX * 


1 1»LABLCYY»2 » LABLPTH* 1 ) 

30 CALL DDIPLT ( 1 » ING,NCYY,XOAtA»SAVPTH,CYMIN,CYMAX,PTHMIN*PTHMAX* 
1 1, LABLCYY, 2, LABLPTH»0> 

IF (NCYY «LT«2) GO TO 40 

CALL DDIPLT < 0 , ING, ENDCYY , ENDPE »CYMIN»CYMAX* EM IN* EM AX* 

1 l,LABLCYY,8,LA6LEN, i> 

CALL DDIPLT (0* ING,2,ENDCyY,ENDIE,CYMIN,CYMAX,EMIN»EMAX* 


1 1,LABLCYY, 
CALL DDIPLT (Of 
1 I , LABLCYY * 
CALL DDIPLT ( 0 » 
1 1, LABLCYY, 

40 CALL DDIPLT (0* 
1 1, LABLCYY* 

CALL DDIPLT ( 0 * 
l 1* LABLCYY, 
CALL DDIPLT (0, 
1 1, LABLCYY, 


8, LABLEN, 2) 

ING,2»ENDCYY,ENDKE,CYMIN,CYMAX*EMIN*EMAX, 

8 , LABLEN » 3 ) 

tNG,2»ENDCYY*ENDTE,CYMlN,CYMAX , EM IN, EM AX, 

8, LABLEN ,4) 

ING,NCYY,XDAtA,SAVPE,CYMIN,CYMAX,EMIN»EMAX» 

8 , LABLEN , 0 ) 

ING, NCYY, XDA T A,5AVIE,CYMIN,CYMAX, EMIN, EMAX, 
8 , LABLEN , 0 ) 

ING,NCYY,XDATA,SAV<E,CYMlNf CYM AX, EM IN, EM AX, 
8, LABLEN, 0) 


CALL DDIPLT (1, ING, NCYY, XDAtA,SAVTE, CYM IN, CYM AX, EM IN, EMAX, 
1 1, LABLCYY, 6, LABLEN, 0) 

RETURN 

END 


) 

K. 





APPENDIX G 


Comparison of Computer Plots Produced by the Polar Coordinate 
Rotating Gas Simulator of Appendix F and an Earlier Rectangular 
Coordinate Code. 

The initial conditions for the two codes are identical. The 
galaxy’s mass consists of a 5% dynamic gaseous component of initial 
maximum radius of 10 cells and a 95% constant stellar component which 
is a uniform density sphere of radius 6 cells. The stellar component 
is represented by the application of an analytically computed constant 
central force. The gaseous component has an initial density distribution 
which varies as [1 - (r/^) 2 ] 1 ^ 2 , where r and t q are the radius and 
maximum radius respectively. The gaseous component has an initial cellular 
specific internal energy distribution which corresponds to a constant times 
the minimum velocity dispersion required to satisfy the Toomre stability 
criterion for a particulate system. Initially, the cells ?iave zero radial 
velocity and radially balancing angular velocities. Both runs were made 
on a 32 x 32 mesh (or polar equivalent) while later runs will be made on a 
64 x 64 mesh. 

The superiority of the polar coordinate code is most readily apparent 
from the last two plots of each run, which show the variation of total 
angular momentum and energy with time. Although the rectangular coordinate 
code conserves total energy, it rapidly looses angular momentum. This loss 
of angular momentum is accompanied by a non-physical decrease in kinetic 
energy and a non-physical increase in internal energy of equal magnitude. 

The later causes a non-physical increase in pressure which along with some 
rectangular grid effects drives the particles rapidly outward. (Because of 
the rapid non-physical rise in internal energy, the scales of the temperature 

^Toomre, Alar: On the Gravitational Stability of a Disk of Stars. 

Astrophys. J., vol. 139, no. 4, May 15, 1964, pp. 1217-1238. 


and pressure plots from the rectangular coordinate code are 11 times smaller 
than those scales for the polar coordinate code). The last two plots of the 
polar code show that it conserves angular momentum exactly and that for the 
first 60 cycles (about one half of one rotation) no heating is evident. 

(The printouts show a slight heating.) The slight outward drift of particles 
near the maximum radius in the polar code is due partly to the physical effects 
of the pressure gradient and partly to a slight instability in the radial 
acceleration, which is currently being investigated. 
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